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SUMMARY 


The main objective of this investigation is to develop a two-equation 
turbulence model for dilute vaporizing sprays or in general for dispersed two- 
phase flows including the effects of phase changes. The model that accounts 
for the interaction between the two phases is based on rigorously derived 
equations for the turbulence kinetic energy (K) and its dissipation 
rate (e ) of the carrier phase using the momentum equation of that phase. 
Closure is achieved by modeling the turbulent correlations, up to third order, 
in the equations of the mean motion, concentration of the vapor in the carrier 
phase, and the kinetic energy of turbulence and its dissipation rate for the 
carrier phase. The governing equations are presented in both the exact and 
the modeled forms. 

It is assumed that no droplet coalescence or breakup occurs. This 
implies that the droplets are sufficiently dispersed so that droplet 
collisions are infrequent. The droplets are considered as a continuous phase 
interpenetrating and interacting with the gas phase, and are classified into 
finite-size groups. Further, constant properties for both the carrier fluid 
and droplets are assumed. 

The Eulerian approach adopted here leads to two sets of transport 
equations, one set for the carrier phase (primary air issuing from the pipe 
plus the evaporated material) and the other for the droplets. These equations 
are coupled primarily by three mechanisms, the mass exchange, the displacement 
of the carrier phase by the volume occupied by droplets, and the momentum 
interchange between droplets and the carrier phase. 

An expression for calculating the turbulent Schmidt number of the 
droplets (the ratio of droplet diffusivity to fluid point diffusivity) is 
developed via comparison with the experimental data (Snyder and Luraley, 1971, 
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and Wells and Stock, 1983). 

The governing equations are solved numerically using a finite-difference 
procedure to test the presented model for the flow of a turbulent axisymmetri 
gaseous jet laden with either evaporating liquid droplets or solid 
particles. The predictions include the distribution of the mean velocity, 
volume fractions of the different phases, concentration of the evaporated 
material in the carrier phase, turbulence intensity and shear stress of the 
carrier phase, droplet diameter distribution, and the jet spreading rate. 
Predictions obtained with the proposed model are compared with the data of 
Shearer et al . (1979) and with the recent experimental data of Solomon et al . 
(1984) for Freon-11 vaporizing sprays. Also, the predictions are compared 
with the data of Modarress et al . (1984) for an air jet laden with solid 
particles. The predictions are in good agreement with the experimental data. 
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1 . 0 INTRODUCTION 


1.1 The Problem Considered 

Dispersed flow is a particular class of two-phase flows, characterized by 
the dispersion of solid particles, liquid droplets, or gas bubbles in a 
continuous fluid phase. Different flow regimes may be encountered. Of 
particular interest here is the case where liquid droplets occupy a small 
fraction, less than 1%, of the total volume of a gas-droplet mixture. This 
spray regime (Fig. 1-1), which has been termed "thin spray" (O'Rourke, 1981) 
or "dilute spray," (Mostafa and Elghobashi, 1984) is important in a variety of 
applications. Steam generators, nuclear reactors, cooling systems, premixed- 
prevaporized gas turbine combustors, diesel-engine sprays, spray-cooling and 
spray-drying systems, and rocket plumes are some examples. Understanding the 
interaction between the particles and surrounding gas is essential for 
predicting dispersed two-phase flows. 

A quantitative definition of "dilutness" in turbulent two-phase flows is 

not readily available. For laminar flow the dilutness requires that the 

center-to-center distance between particles should be larger than 2(a+6_) 

b 

where a is the particle radius and <$, is the thickness of the boundary layer 

b 

around that particle. The experimental data of Tsuji et al. (1982) indicates 
that the fluid dynamic force on a suspended particle can be assumed to be the 
same as that on a single particle if the interparticle spacing is not less 
than three particle diameters. This restriction gives an upper limit for the 
volume fraction of the particles of 1% to satisfy the diluteness assumption. 
For turbulent flow, another parameter plays a more significant role. This is 
the ratio of the aerodynamic response time to the time between collisions 
(Crowe, 1981). This ratio depends on the particles' loading ratio, the 
relative velocity between particles and gas, and the gas velocity gradient. 
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Figure M THE FLOW CONSIDERED 
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If this ratio is less than unity, a particle has time to respond to the local 
gas velocity field before the next collision so its motion is dominated by 
aerodynamic forces and the particle collisions can be neglected. Using this 
restriction the particle number density, or the volume fraction, for a dilute 
spray can be calculated. 

In the dilute spray regime the interactions between droplets is 
neglected. This implies that the droplet coalescence or break-up does not 
occur. The droplets may exchange mass, momentum, and energy with the gas and 
for dilute spray the exchange functions for isolated droplets can be used. 

The droplets are classified into finite-size groups and each group is 
considered as a distinctive phase. Further, constant properties are assumed 
for all the phases to avoid the density fluctuations of the carrier phase at 
this stage. 


1.2 Previous Work 

1.2.1 Fundamentals of Two-Phase Flow 

The simplest analytical approach for calculating the properties of dilute 

suspensions of two-phase flow is to assume dynamic equilibrium, where the 

particles and gas velocities are equal at each point in the flow. The 

suspensions can then be considered as a single homogeneous fluid that is 

treated exactly as a single-phase flow. The mixture properties are those 

based on the continuum mechanics that apply to molecular mixtures. The 

infinitely fast interphase transport between the phases is the basic premise 

of that approach. The equilibrium assumption is valid for small values of 

Stokes number, less than 10 - ^, and small values of particle/fluid material 

o 

density ratio, less than 10 , (DiGiacinto et al., 1982). Stokes number is the 
ratio of the particle relaxation time to the characteristic time of the 
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surrounding fluid. These two restrictions are not satisfied for the flow of 
gaseous phase laden with liquid droplets or solid particles. Accordingly, the 
local equilibrium approximation leads to unrealistic results for that type of 
flow (Shuen et al. , 1983). In spite of the inaccuracies in that approach, it 
has been used by some workers (Shearer et al., 1979; Michaelides, 1984; and 
Kamimoto and Matsuoka, 1977). 

On the other hand, if both the density ratio and the Stokes number are 
large, the particles will not be able to respond to the changes in the carrier 
phase. In this case the difference between the velocities of the phases can 
not be neglected and each phase should be treated separately (Crowe, 1982). 

There are two approaches to handle the carrier phase in the separated 
flow models, depending on the mass loading ratio, which is defined as the mass 
flow rate of the particles to that of gas. If this ratio is small, less than 
0.1, the velocity field of the carrier phase is not affected by the presence 
of the particles while the motion of the particles is determined by the gas 
flow properties (Rudinger, 1965). In this case the governing equations of the 
carrier phase have no extra terms but rather they are identical to the well- 
known Navier-Stokes equations. This approach is referred to as one-way 
coupling (DiGiacinto et al., 1982) from gas to particles only and has been 
used by many workers (Cox and Mason, 1971; Batchelor, 1974; and Boyson and 
Swithenbank, 1979). On the other hand, if the mass loading ratio is high, the 
particles may modify the gas flow field significantly. In this case the 
particles are regarded as source of mass and momentum for the carrier phase. 
This approach is referred to as two-way coupling (DiGiacinto et al., 1982; and 
Crowe, 1982) from fluid to particles and vice versa. 
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There are two main approaches to handle the dispersed phase in the 
separate two-phase flow models, namely the Lagrangian and the Eulerian 
approaches. In the Lagrangian approach the dispersed phase is treated by 
solving Lagrangian equations of motion for the particles with a prescribed set 
of initial conditions. Once the flow properties of the particles are known, 
the interface quantities between the two phases can be calculated. In the 
Eulerian approach the dispersed phase is treated as an interacting and 
interpenetrating continuum. In that approach the governing equations for the 
two phases are quite similar to the well-known Navier-Stokes equations. These 
equations are coupled primarily by three mechanisms, the mass exchange, the 
displacement of the carrier phase by the volume occupied by particles, and 
momentum interchange between particles and the carrier phase. Many two-way 
coupling studies are presented in the literature, based either on the 
Lagrangian or Eulerian approaches (Elghobashi and Megahed, 1981; Yeung, 1982; 
Abbas et al., 1981; and Crowe et al., 1977). 

Of most importance, the continuum assumption must be justified when using 
the Eulerian approach. Batchelor (1974), Lumley (1978a), and Marble (1962) 
have discussed the continuum concept for the dispersed phase. In summary, the 
particles must be sufficiently small in order that a volume element, small 
compared to the Kolmogoroff microlength scale, n, contains a large number of 
particles. Thus a statistical average concerning the behavior of the 
particles can be made within this volume element. This requires that the 
average separation distance between the particles is at least one order of 
magnitude smaller than h. Hinze (1972) stated that the continuum assumption 
has proven to be applicable also to situations that do not strictly meet that 
condition. Others (Crowe, 1982; Soo, 1967; and Yeung, 1978) showed that most 
practical physical systems involving gas-particle mixtures satisfy the 
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continuum assumption. We may refer, amongst others, to theoretical 
investigations by Marble (1970), Buckingham and Siekhaus (1981), Pourahmadi 
and Humphrey (1983), Rizk and Elghobashi (1984), and Mostafa and Elghobashi 
(1983) who used the Eulerian approach to study different flow conditions of 
.two-phase flows. Early work based on Lagrangian equations of motion are due 
to • El-Banhawy and Whitelaw (1980), Mongia and Smith (1978), Shuen et al. 
(1983), El-Kotb et al . (1983), and El-Emam and Mansour (1983). 

Arguments over the advantages and the disadvantages of Eulerian and 
Lagrangian approaches persist in the literature. The Eulerian approach models 
can easily incorporate particle diffusion effects since the randomness of the 
particulate phase is accounted for by the way of the formulation. This 
approach can be extended easily to multidimensional flows. Numerical 
instabilities, false diffusion, and large storage requirements are the most 
serious disadvantages of that approach. However, the use of advanced digital 
computers and the ability to overcome the numerical problems (for example, by 
choosing a suitable higher order finite-difference scheme) alleviate most of 
these disadvantages. The Lagrangian approach exhibits no numerical diffusion 
but the particle dispersion must be incorporated through an empirical 
diffusion velocity or more expensive Monte Carlo methods (Chen and Crowe, 
1984). Durst et al. (1984) showed that the Lagrangian approach, in cases 
where the particle loadings are high, is inferior to the Eulerian approach. 

The Lagrangian approach calculations require interpolation between the meshes 
since gas and particle properties are strongly coupled. In any case it 
requires a toilsome computation for the source terms. If the interpolation 
process is too crude, Aggarwal et al. (1983) have shown that errors of the 
same order as the diffusion error in the Eulerian approach will be 
encountered. Sirignano (1983) argued that the droplet properties should not 
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be averaged over the numerical cell as suggested by Dukowicz (1980), but 
rather a linear interpolation should be made. In the present work attention 
will be restricted to a formulation following the Eulerian approach. 

The previous fundamental studies of the various aspects of two-phase flow 
are concentrated on either the effects of various factors on the flow around a 
single particle or on the governing equations of the dispersed phase. Fuchs 
(1964) and Torbin and Gauvin (1959, 1960, 1961) did an extensive survey about 
the dynamics of single particles. Those studies are very important to 
fundamentally understand the two-phase flow. Most of the recent publications 
in this regard will be discussed in section four. 

On the other hand, several phenomenological attempts have appeared in the 
literature to derive equations governing the macroscopic behavior of dispersed 
two-phase flow. The equations cited most frequently are those of Drew (1971), 
Kalinin (1970), Whitaker (1973), Gray (1975), Panton (1968), and Soo (1967). 
Other deriviations include those by Nigmatulin (1967), Owen (1969), Rietema 
and Van Den Akker (1983), Buevich and Markov (1973), and Jackson and Davidson 
(1983). The resulting equations differ in various ways such as the 
formulation for the pressure gradient term, the nature of the momentum, source 
terms, or the proper coupling between the two fields. 

Buevich and Markov (1973) obtained the conservation of mass, momentum and 
moment of momentum for the two interpenetrating and interacting continua. All 
the unknowns in the governing equations are expressed in terms of mean 
stresses acting at the surface of an individual suspended sphere. Crowe 
(1980) used the control volume, or Reynolds transport theorem approach, to 
derive the continuity and momentum equations for a flowing vapor with 
suspended burning, evaporating, or condensing droplets. 
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Solbrig and Hughes (1975) derived the momentum equations and mechanical 
constitutive equations that are required to describe transient, two-phase, 
single-component evaporating and condensing flows. Momentum field balance 
equations were derived for each phase on the basis of a seriated-continuum 
approach. A seriated continuum is distinguished from an interpenetrating 
medium by the representation of interphase friction with velocity differe nces 
in the former and velocity gradients in the latter. A two-phase mixture is an 
example of a seriated continuum, whereas a mixture of gases is an example of 
an interpenetrating continuum. The seriated continuum also considers embedded 
stationary solid surfaces such as that which occurs in nuclear reactor 
cores. There are some undetermined numerical coefficients that appeared in 
the momentum equations of Solbrig and Hughes (1975). These coefficients must 
be determined for the different flow regimes and geometry. 

Panton (1968) formulated the flow properties for the non-equilibrium two- 
phase flow of a gas-particle mixture. The conservation equations of continuum 
fluid mechanics are assumed to apply to the flow field locally, both within 
the particles and through the gas. Control volumes for each phase are defined 
and integral forms of the conservation equation are applied. By inspecting 
the equations, the proper area-averaged properties are defined so that they 
are meaningful terras in the physical conservation laws. Because the detailed 
flow is inherently unsteady, it was necessary to take the time average of the 
equations. Thus, the dependent variables of the the final conservation 
equations were area-time-averaged properties. New terms, even in laminar 
flow, appeared in the momentum equations and were called the area-averaged 
Reynolds stresses. The Reynolds stresses attributed to the fluctuations in 
the gas velocity occur because of the presence of the particles. Every time a 
particle passes the point under observation, a fluctuation in the gas velocity 
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occurs . Panton commented that these unknowns are the price to be paid for the 
details of the flow. 

Delhay (1980) surveyed two-phase flow modeling. He discussed the 
different types of averaging — time, area and volume — in two-phase flow. 
Also, the different two-phase flow and single fluid modeling were reviewed. 
Delhay concluded that more work is needed to build the bridge between the 
Lagrangian behavior of a particle and the Eulerian form of the constitutive 
terms entering the averaged balance equations. 

Drew (1971) derived averaged field equations for two-phase media. He 
treated the separated surfaces between the two-phase media as transition 
regions where the material properties have jump discontinuities. Postulating 
the laws of balance of mass, linear momentum, angular momentum, energy, and an 
entropy inequality, jump condition laws for each phase were derived. Solving 
the differential equations and jump conditions, exact expressions for the 
field quantities involved were found. Drew also defined and related the 
appropriate average variables for each phase involved. He commented that for 
any particular problem, his averaged field equations must be supplemented by 
constitutive equations, which is not a simple task. Ishii (1975) discussed 
the way of averaging used by Drew (1971), specifically the two integrals over 
both space and time domains. Ishii commented that it is not quite convincing 
why these four integrations are necessary to develop meaningful macroscopic 
field equations. He pointed out that the time and space differential 
operators in the averaged fields represent finite difference operators in the 
physical interpretations. 

Ishii (1975) presented a detailed discussion on the formulation of 
various mathematical models of two-phase flows based on the conservation laws 
of mass, momentum and energy. He considered the local instant formulation and 
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the time-averaged macroscopic models. He presented the two-fluid model, which 
is formulated by considering each phase separately. Thus, the model is 
expressed by two sets of conservation equations of mass, momentum, and energy 
with interaction terms appearing in the field equations. His formulation has 
the advantage of treating large and small particles alike, with averaging 
carried out across the interface. Ishii's formulation simplifies the 
treatment of the dispersed phase by introducing a duality of discrete nature 
and distributive representation. The discreteness is accounted for via 
treating the virtual mass and unsteadiness of flow field of each finite size 
particle. The distributive nature of the particle cloud is accounted for by 
taking an elementary volume consisting of a sufficiently large number of 
particles. Ishii also considered the diffusion model, which is formulated by 
considering the mixture as a whole. Thus, it is expressed in terms of three 
mixture conservation equations of mass, momentum, and energy with one 
additional diffusion equation. 

Sha and Soo (1978) discussed the basic concepts for the rigorous 
formulation of a system of a single-component fluid in two phases. They 
pointed out that the direct extension of continuum mechanics is inadequate 
because of the mutually exclusive nature of the phases in a multiphase 
system. Multiphase mechanics have their own distinct regime with additional 
inertial and viscous interaction terms, applied to mixtures of phases that are 
separated by interfaces and are mutually exclusive. This is in contrast to 
the field equations of mixtures based on continuum mechanics, which directly 
apply to molecular mixtures where the phases coexist at the same points in 
space. Boure (1979), Crowe (1978), and No (1982) argued that the equations of 
Sha and Soo (1978) are inconsistent and not valid even in one-dimensional 
situations . 
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In section two, the governing equations of dispersed two-phase flow are 
presented and compared with other equations in the literature. 

1.2.2 Turbulent Evaporating Sprays 

Modeling of evaporating and combusting sprays is an extremely difficult 
problem due to the complex physical and chemical phenomena encountered in this 
type of two-phase flow. A substantial number of reviews of this problem have 
appeared in the literature. The recent reviews of Law (1982), Faeth (1977, 
1983), Labowski and Rosner (1973), and Sirignano (1983) discussed the previous 
work on the different phenomena associated with the spray evaporation and 
combustion problem. The present study will be restricted to evaporating or 
nonevaporating dilute sprays. 

Krestein (1983) has analyzed a simple model of an evaporating spray to 
predict the probability density function (pdf) of vapor concentration within 
the spray. The model assumes that the droplets deposit linear streaks of 
vapor as they traverse the motionless host gas, and that the vapor diffuses 
radially from these streaks. Since it neglects droplet collisions, 
saturation, and related effects, the model is applicable primarily to dilute 
sprays. The results of this analysis can be used to estimate droplet 
vaporization rates from experimentally measured pdf's of concentration. 
Therefore the individual-droplet processes could be linked to fluctuating 
ambient conditions in spray simulation codes. 

O'Rourke and Bracco (1980) developed a numerical model for turbulent 
dense sprays. The model is two-dimensional unsteady and uses atomization 
experimental results as nozzle exit boundary conditions and a stochastic 
algorithm to compute droplet events, including collisions and coalescence. 
Westbrook (1976) presented a numerical solution technique for the spray 
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equation for' a type of stratified charge internal combustion chamber. He 
neglected the gas entrainment by the droplets and adopted the dilute spray 
approximations. The gas motion was assumed to be consisted of a rotational 
swirl with a constant angular velocity. Axial and radial components of the 
gas velocity were assumed to be identically zero. 

Martin, elli et al. (1983) used O'Rourke's model after considering a 
K-ie submodel for gas turbulence to predict the data of Wu et al. (1984). 
Agreement is good with mean quantities but the computed standard deviation of 
the drop velocity distribution is generally smaller than the measured one. 
Although the effects of turbulence on the droplet motion is considered in the 
model, the direct effects of the droplets on the gas motion are neglected. 

Yeul et al. (1982) have since reported measurements in evaporating 
kerosene sprays from a twin-fluid injector in a co-flowing stream. 
Measurements of droplet size were undertaken using a laser tomographic light- 
scattering technique while mean velocities were measured using LDA. They did 
not measure the turbulence characteristics or the droplet/velocity 
correlations which are needed for the theoretical models evaluation. Wu et 
al. (1984) reported LDV measurements for the distribution function of the 
axial and radial components of the droplet velocity at various radial and 
axial locations within steady sprays under the conditions of direct fuel 
injection in internal combustion engines, but at room temperature. The 
measurements were taken within 300 to 800 nozzle diameters from the nozzle 
exit. 


1.2.3 Turbulence Mathematical Models 

Computational models are a very useful tool for a better understanding of 
the features of the two-phase flow, considering the inability of the 
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analytical methods and the difficulty of experimental investigations. 

Vasillev (1969) reviewed the development of the two-phase flow, relying 
chiefly on the research that has been done in the Soviet Union. He cited the 
papers on the governing equations for laminar and turbulent flows as well as 
those on the effects of the dispersed phase on the turbulence intensity and 
the spectrum of turbulence. He concluded that the presence of small suspended 
particles leads to more rapid damping of the turbulent energy under isotropic 
flow conditions. In the case of large values of density ratio it also causes 
a noticeable distortion of the turbulence energy spectrum and a decrease of 
its micro scales in comparison with the case of single-phase flow. 

Rakhmatulin (1956) (cited by Vasiliev, 1969) suggested that the motion of 
the mixture should be treated as an interpenetrating motion of several 
continua. The equations of motion are written separately for each phase, and 
the interaction between the phases is taken into account by considering the 
interaction forces that appear as internal forces for the whole system. The 
governing equations of Rakhmatulin (1956) were used by Bondarenko and 
Shaposhnikova (1980) to analyze flow regimes in channels of different 
shapes. Those equations were also used by Vasil'kov (1976) to predict a 
turbulent submerged jet containing an admixture of solid particles. 

Michaelides (1984) analyzed the gas-solid two-phase pipe flows using the 
mixing length hypothesis. The mixture was taken to be a homogeneous fluid of 
variable density across the pipe cross-section. These two assumptions make 
the solution very restricted to dilute suspensions with comparable densities 
between the solid and the gas. 

Buckingham and Siekhaus (1981) described a K-e turbulence model that 
allows for effects of particles on turbulence properties. The model is 
applied to flows containing small solid particles, considering added 
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dissipation due to particle interactions with the carrier phase in the 
governing. Buckingham and Siekhaus did not compare the performance of their 
model with any experimental data. The predictions suggest a damping of the 
turbulence motions primarily because of inertial effects. Nagarajan and 
Murgatroyd (1971) presented an analytical model for turbulent two-phase fully- 
developed pipe flow. They assumed linear shear stress in the radial dir ect i° n 
and introduced several phenomenological coefficients in the model. This made 
their model inapplicable to other two-phase flow problems. 

Kramer and Depew (1972) developed a one-dimensional model for a fully 
developed two-phase turbulent pipe flow. In their solution they expressed the 
velocity fields in terms of various empirical coefficients and assumed a 
linear mixing length to express the turbulent correlations. This has again 
made the application of their model to any other problem very difficult. 

Genchev and Karpuzov (1980) have proposed a turbulence model for fluid-~ 

particle flows in which the effects of particles on the turbulence transport 

equations are considered. They assumed that the mean velocity of the 

particles is equal to the fluid mean velocity and neglected the fluid-particle 

turbulent correlations existing in the time-averaged equations (Elghobashi r -and 

Abou-Arab, 1983). Genchev and Karpuzov predicted a fully developed pipe flow 

laden with solid particles of density ratio of order 10 and volume 

—3 

concentration of order of 10 . They did not compare their predictions with 

experimental data to evaluate the capabilities and the limitations of their 
model . 

Danon et al. (1977) described a K-L model for two-phase jets. The length 
scale (L) was not modified from the value appropriate for a constant density 
single phase jet; however, a term representing the added dissipation due to 
the presence of particles was included in the governing equation for K. The 
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model was evaluated using the data of Hetsroni and Sokolov (1971) for a round 
jet containing oil droplets. The predictions using the basic model were not 
in good agreement with these measurements. The comparison between the 
predictions and measurements was improved by multiplying the rates of 
production and dissipation of K by a coefficient that was a strong function of 
the void fraction. These authors commented that there is a substantial and 
unexplained influence of particles on turbulence properties of jets, even at 
low particle concentrations. 

Melville and Bray (1979) described a model for particulate flow, with 
small interphase slip, employing constant eddy dif fusivities for momentum and 
particle transport. The predictions were evaluated using the measurements of 
Laats and Frishman (1970a and 1970b) for a round jet containing powders of 
various sizes. 

Pourahmadi and Humphrey (1983) proposed a mathematical model for dilute 
suspensions of two-phase flow based on the single-phase K-g model. They 
considered the direct effects of the particle's sharing the same control 
volume with the gas on the governing equations of K and g . These authors 
neglected all the third-order correlations without justification and used 
stokes drag coefficient although the particle Reynolds number is generally 
greater than unity in two-phase flow. They also used Peskin's formula for the 
calculations of the particle's Schmidt number. This formula was tested by 
Elghobashi et al. (1984) for glass particles and produced unrealistic values 
for the particle's Schmidt number, negative or zero. 

Elghobashi and Abou-Arab (1983) proposed a two-equation turbulence model 
for incompressible dilute two-phase flow which undergoes no phase changes. 
Using this model, Elghobashi et al. (1984) predicted the turbulent 
axisymmetric gaseous jet laden with uniform size solid particles. They 
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achieved good agreement with the experimental data of Modarress et al. 

(1984). This model has been extended by Mostafa and Elghobashi (1985a) to 
include the effects of phase changes. 

Shuen et al. (1983) evaluated the performance of the available Lagrangian 
methods for predicting the dispersed phase behavior by comparing the results 
with measurements of particle-laden jets. They considered only dilute 
suspensions of solid particles, and hence concluded that the effects of the 
particles on the turbulence quantities are almost negligible. This allowed 
them to recommend the use of the conventional K-e model for two-phase flows 
without any modifications. Shuen et al. (1983) indicated that the suggested 
method by Gosman and Ioannides (1981) for calculating particle trajectories, 
the "stochastic or Monte-Carlo method," in contrast to other methods, provides 
good predictions over their data base. In this method the isotropic turbulent 
gas velocity field is split into mean and fluctuation. The mean value is 
obtained from the solution of the mean equations while the fluctuating one is 
estimated from random sampling of a Gaussian distribution of the kinetic 
energy of turbulence. The Monte-Carlo method requires selection of 
characteristic eddy length and time scales. Shuen et al. (1983), following 
Gosman and Ioannides (1981), assumed that the eddies are uniform and their 
size is proportional to a turbulent length scale, £ e , given by 

i = C K 3/2 /e 1.1 

e e 

Shuen et al. (1983) used the value of 0.16 for the constant C e , which was 
suggested by Gosman and Ioannides (1981) who later changed this value to 0.31 
to get a better agreement with the experimental data (Crowe, 1982). Crowe, 
(1982) recently argued that the value of C g should be 0.46 to give a good 
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agreement with the experimental results. This is unfortunate since Shuen et 
al. (1983); Gosman and Ioannides (1981); and Crowe, (1982) have used the same 
experimental data of Snyder and Lumley (1971) to obtain the value of the 
empirical coefficient C e . 

In conclusion, after the literature on two-phase turbulent flow was 
examined, it was found that there is still no complete mathematical model of 
this class of flows comparable with the model of single-phase turbulent 
flows. The main objective of this study is to develop such a model. 

1.2.4 Turbulent Two-Phase Jet Flows 

A turbulent nonreacting gaseous jet laden with solid particles or 
evaporating droplets is a relatively simple flow that allows the study of the 
interactions between the two phases and the turbulent dispersion of the 
discrete phase. Most of the previous measurements (Rajani, 1972; Hetsroni and 
Sokolov, 1971; and Laats and Frishman, 1970a) of the structure of particle- 
laden jets considered the effects of the dispersed phase on the continuous 
phase properties. Abramovich (1970) and Goldschmidt and Eskinazi (1966), 
discussed the effects of the dispersed phase on the structure of a turbulent 
gaseous jet. They showed that the particle concentration profiles in a two- 
phase jet are narrower than the gas velocity profiles. This behavior was 
explained later by Elghobashi et al. (1984). 

Levy and Lockwood (1981) and Laats and Frishman (1970a) found that the 
gas mean velocity profiles in a two-phase jet are narrower than those of the 
clear jet. Modarress et al. (1984) and Girshovich et al. (1981) further 
showed that the solid or liquid particle velocity is higher than the gas 
velocity in the developed region of the jet. Other studies (Al-Taweel and 
Laundau, 1977; and Laats and Frishman, 1973) showed that the turbulent energy 
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level decreases with the increase of the suspension of particles into a jet. 
The effect of small droplets of cottonseed of 13 pm average diameter on the 
flow structure of an axially symmetrical turbulent air jet has been studied by 
Hetsroni and Sokolov (1971). They found that in the two-phase jet, the 
velocity spread and the turbulence intensities were reduced in comparison with 
the single-phase jet. They also found that even at low volumetric droplet 
loadings, the jet was narrower than single-phase air jet. At a high loading, 
the jet spread was wider upstream and narrower in the downstream region. The 
intensity of velocity fluctuations was reduced throughout the jet. Hetsroni 
and Sokolov (1971) measured time-averaged and fluctuating longitudinal 
velocities by means of a hot-wire anemometer. The probe was not calibrated in 
a two-phase flow but the authors stated that the calibration curves obtained 
in single-phase flow could be used for two-phase flow with minor 
corrections. Rajani (1972) pointed out the uncertainties regarding the 
calibration of probes used in dust-laden flows that may lead to an 
overestimation of the measured quantities. Therefore, their results should be 
viewed with caution as pointed out by Melville and Bray (1979). 

Field (1963) and Subramanian and Ganesh (1982a and 1982b) have provided 
evidence of the overall effect of solid particles on a dust-laden jet by 
measuring the rate of ambient air entrainment by the jet. Field considered 
lycopodium powder of 30 pm size while Subramanan and Ganesh used sand 
particles of uniform size of 150-180 pm. They found that the entrainment was 
affected by particle size, density and mass loading ratios, and the dispersed 
phase initial conditions. Subramanian and Ganesh showed that the presence of 
particles increases the entrainment rate. Since their measurements are in the 
developing region (z/D^_ 7), where z/D is the axial distance to nozzle 
diameter ratio, and no nozzle exit conditions are reported, it is difficult to 
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analyze that data. 

Hayashl and Branch (1980) measured the concentration profiles of 
particles in axisymmetric jets. The measurements were made by seeding 1%, 3%, 
and 5% by weight of 24 pn mean diameter spherical flash ash particles into 
jets at Mach numbers of 0.2, 0.8, and 1.0. The particle concentration 
profiles showed that particles concentrate on the axis of the jet at the exit 
of the nozzle and the profiles are highly influenced by the initial 
conditions. This observation is expected since the measurements were done in 
the developing region, which is highly affected by the nozzle exit conditions. 

Zuev and Lepeshinskii (1981) studied the two-dimensional steady isobaric 
two-phase jets. They considered the effects of particle-particle interaction 
on the governing equations from an analogy with the kinetic theory of gases. 
They adopted the mixing length hypothesis to close the set of equations. Zuev 
and Lepeshinskii did not compare the predictions using their model with any 
experimental data to test their approach. Vasil’kov (1976) added terms that 
take into account the interaction of the phases to the governing equations of 
single-phase gas dynamics to predict a turbulent submerged jet containing an 
admixture of solid particles. He assumed that the radial velocities of the 
particles are equal to those of the gas and adopted the mixing length 
hypothesis to close the set of the governing equations. Vasil 'kov obtained a 
reasonable agreement between the predictions using his phenomenological model 
and the data of Laats and Frishman (1970a and 1970b). 

Popper etal. (1974) studied the motion of oil droplets of 50 pm, in a 
round air jet using LDV. It was found that at the jet axis, the droplet 
velocities are 5-9% lower than the corresponding velocities in a single-phase 
air jet. In the developed region the droplet velocities were higher than the 
air velocity at the same location. 


19 


The effects of spherical glass particles of 110 un average diameter on 
the flow structure of an axially symmetrical turbulent air jet has been 
studied by Rajani (1972) for various solid/air loading ratios, from 0 to 1.3 
kg/kg. Time-averaged velocity measurements were performed by laser Doppler 
anemometer (LDA) and the particle concentration measurements by a scattered 
light technique. Rajani devoted a large part of his work to the development 
of the experimental techniques, their accuracy, and limitations. Therefore, 
he reported very limited data, especially for the dispersed phase. 

Yuu et al. (1978) examined the distribution of _ concentrations of the dust 
particles of the average diameter on a mass basis of 15 and 20 |n In a round 
jet. The flow was highly dilute, since solid volume fraction in the injected 
flow was in the range of 0.4 - 2xl0 - ^. The measurement of particle 
concentrations was performed with a photoelectric dust counter and the mean 
velocity was measured with a pitot-static probe. Using the concentration 
data, they indicated that the particle diffusivity decreases with the increase 
of the particle inertia and in general it is smaller than that of fluid scalar 
quantities. Goldschmidt and Eskinazi (1966) measured the concentration of the 
liquid droplets of the average diameter on a weight basis of 3.3 un in a two- 
dimensional jet. They indicated that the droplets' mass tends not to diffuse 
more than the fluid momentum. 

Girshovich et al . (1981) and Laats and Frishman (1973) investigated 
experimentally the effect of solid particles on both the mean and turbulence 
axial velocity components of an air jet using LDV techniques. They studied 
the effect of the circular tube diameter, initial velocities, initial 
particles/air mass loading ratio, and particle diameters on the jet 
performance. Laats and his coworkers found less rapid decay of centerline 
velocity and a reduced velocity spread of the jet with the increase of the 
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solid loading. They also found that the increase in the mean axial velocity 
and the decrease in the turbulence intensity at the centerline for the carrier 
phase depend on the loading ratio and the diameter of the particles. 

Unfortunately Laats and his coworkers did not measure the initial conditions 

at the nozzle exit and did not report the material density of the particles. 
This fact renders their data inadequate for the evaluation of turbulence 
models . 

Levy and Lockwood (1981) measured fluid and solid phase mean and 
fluctuating velocities in a. round gaseous jet using LDV techniques. They 
studied sand particles ranging in size from 215 to 1060 pm with sand to air 
mass ratios ranging from 1.14 to 3.5. Levy and Lockwood found that, relative 
to the pure gas phase, the axial turbulence intensity was reduced by 
introducing particles in the size range of 180-500 pm and was increased when 
the particles are in the range of 500-1200 Mn. But again they did not report 
the nozzle exit conditions. Modarress et al. (1983) reported much needed 
experimental data to help understanding the behavior of two-phase turbulent 
jets and to validate the theoretical models for these flows (Fig. 1-1). They 
investigated the effects of 50 Wn and 200 pn glass beads on the mean air 

velocity and found that the turbulent stresses for mass loading ratio varies 

from 0.32 to 0.85. Modarress et al. found that the increase in the centerline 
mean air velocity and the diminishing of the turbulence quantities are 
proportional to the loading ratio, the particle diameter, and the initial 
conditions of each phase. 

Shearer et al. (1979) measured the mean velocity, velocity fluctuations, 
and Reynolds stress of single-phase constant density jets, as well as those of 
an evaporating spray (Freon-11). Since their measurements were in the region 
far downstream from the nozzle (170 <_ z/D <_ 510), the flow was highly dilute 
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Figure 1-2 THE FLOW CONSIDERED (SOLID 

PARTICLES OR METHANOL LIQUID 
DROPLETS) 
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and the effects of the dispersed phase on the gas properties were very 
slight. Solomon et al. (1984) measured the flow properties of the carrier 

phase as well as those of the droplets in a turbulent round jet laden with 

Freon-11 spray. They considered the dilute portion of the spray (50 < z/D < 

510) injected into a still air environment in order to provide data useful for 
the evaluation of spray models. They measured all the radial profiles of the 
main dependent variables at 50 nozzle diameters from the exit plane for two 
mass loading ratios of 7.71 and 15.78. This information is essential for 
accurately predicting such flow. 

The last three experiments (Modarress et al., 1983; Shearer et al., 1979; 
and Solomon et al., 1984) are used in the this work to test the proposed 
turbulence model. 

1.3 Summary of Approach 

The present contribution focuses on developing and testing a two-equation 
turbulence model for predicting isothermal steady two-phase flows including 
the effects of phase changes. A set of equations describes the conservation 
of mass, momentum of each phase, vapor concentration, and kinetic energy of 
turbulence and its dissipation rate for the carrier fluid. Closure of the 
time-mean equations is achieved by modeling the existing turbulent 
correlations up to the third order. The model considers turbulent non- 
reacting axisymmetric jet flows laden with evaporating droplets or solid 
particles. This flow regime is a relatively simple flow that allows the study 
of the interactions between the droplets and the carrier phase, the turbulent 
dispersion of the droplets. The radial profiles of the main dependent 
variables are easy to measure in this type of flow, thus it is convenient for 
the turbulence model's validation. 
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In the following sections, the transport equations governing the mean 
quantities are presented first of all, followed by the development of the 
proposed two-equation model. Then the mass and momentum exchange coefficients 
are evaluated and an expression for the particle's Schmidt number is 
developed. The work concludes with an evaluation of the model using the 
recent measurements of Modarress et al. (1983), Solomo n et (1984), and the 
measurements of Shearer et al. (1979). 
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2.0 GOVERNING EQUATIONS OF DILUTE SPRAY 

The purpose of this section is to present the basic equations that govern 
the turbulent dilute vaporizing sprays, and to discuss the problems that their 
solution poses. 

Section 2.1 states the assumptions for this study. Section 2.2 lists 
the time-dependent equations. The discussion then turns to time-averaged 
equations in section 2.3. Finally, the problem of closure is discussed in 
section 2.4. 

2.1 Assumptions 

It is assumed that no droplet coalescensce or breakup occurs. This 
implies that the droplets are sufficiently dispersed so that droplet 
collisions are infrequent. This assumption renders the present study 
restricted to dilute suspensions only. The initial breakup of liquid sprays 
or jets is not considered. It is assumed that the initial profiles of volume 
fractions and velocities are independently specified. Therefore, there is 
two-way coupling between the droplets and the carrier phase. It is assumed 
that the droplets of different sizes constitute different continuous phases. 
This is from the point of view of the "continuum" mechanics of a cloud of 
droplets, apart from the obvious definition of a multiphase system, a mixture 
of phases of liquid droplet and gas (Soo, 1967). Therefore, the continuous 
droplet-size distribution is divided into n intervals; d^ is the average 
diameter for droplets in the k 11 * 1 diameter range. If d s and d^ are the 
smallest and largest droplet diameters, then the sizes are ordered as follows 

d s = d n < d n_1 -— < d 1 = d L 2.1 
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Thus, n different diameter ranges constitute correspondingly n dispersed 
phases and the evaporated mass with the surrounding gas constitute the carrier 
phase. 

It was also assumed that the droplets are sufficiently small in order 
that a volume element, small compared to the Kolmogoroff microlength scale, 
contains such a large number of droplets that a statistical average concerning 
the behavior of the droplets can be made within this volume element. It was 
further assumed that the droplets remain spherical during their entire 
lifetime. This assumption is discussed in detail in section 4.2. Also, it 
was assumed that the mean flow is steady and the material properties of the 
different phases are constant. 

This leads to two sets of transport equations, one set for the droplets 
and the other for the carrier phase, which is defined as the atomizing air 
plus the evaporated material. These equations are coupled primarily by three 
mechanisms, the mass exchange, the displacement of the carrier phase by the 
volume occupied by droplets, and momentum interchange between droplets and the 
carrier phase. The momentum interchange is due to the aerodynamic forces 
exerted on the dispersed phase and the momentum growth resulting from the 
relative velocity between the generated vapor and the surrounding gas. 

2.2 Time Dependent Equations 

As discussed in subsection 1.2.1 many authors derived continuity and 
momentum equations for each phase by performing volume averaging (Sha and Soo, 
1978; Hinze, 1972; and Jackson and Davidson, 1983) or averaging in space and 
time (Panton, 1968; and Drew, 1971). Here the instantaneous, volume-averaged 
equations, in Cartesian tensor notations, are presented based on those of 
(Crowe, 1980; Hinze, 1972; Harlow and Amsden, 1975; and Jackson and Davidson, 
1983). 
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The continuity equation for the carrier phase is 


+ -i* k4k 


2.2 


The term on the right-hand side of Equation 2.2 represents the rate of change 
of the added mass or the source term due to the evaporation process from all 
droplets existing with the carrier phase in the same control volume. This 
term also represents a sink term in the continuity equation of the droplets 
(Equation 2.3). 


The continuity equation for the k^ dispersed phase is 


(P 2 * k ) t + (P 2 * k V k ) . = - m k $ k 


The global continuity is 
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The momentum equations for the carrier phase are 
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The momentum equations for the k dispersed phase are 
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The set of equations (2.2 - 2.6) have (4k + 5) unknowns (3k of droplet 
velocities (V^), k of droplet volume fractions (<j> ) , 3 carrier phase 
velocities (U^), carrier phase volume fraction ($^), and the static pressure 
(P)) and (4k + 5) equations. So it forms a closed set of equations since the 
number of unknowns is equal to the number of equations. 

Using the continuity equations for the different phases (Equations 2.2 
and 2.3), Equations 2.5 and 2.6 as can be rewritten as 
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In the equations above and throughout this work the partial derivatives are 
represented by a subscript consisting of a comma and an index; e.g. , ( ),t 
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The subscripts 1 and 2 denote, respectively, the carrier fluid and dispersed 
phase; the superscript k denotes the dispersed phase; is the velocity 
component of the carrier fluid; is the velocity component of the droplets 
in the k 1 "* 1 diameter range; p and p are the material density and viscosity; P 
is the pressure; $ is the volume fraction; is the gravitational 

acceleration in the i direction; F is the interphase friction coefficient, 

• k 

and m is the evaporation rate per droplet volume. 

Note that the factors 4^ and $ lie outside the pressure gradients in 
Equations 2.5 and 2.6, contrary to what some authors have proposed for 
considering those factors within the gradients (Sha and Soo, 1978). Harlow 
and Amsden (1975), Nigmatulin (1979), Solbrig and Hughes (1975) and Boure 
(1979) argued that ^ and $ should be outside the pressure gradients as in 
Equations 2.5 and 2.6. 

The momentum growth term (m * V^) in Equations 2.5 and 2.6 represents a 
force on the fluid due to the difference between the velocity of vapor leaving 
the droplet surface and that of the carrier fluid. If the flow from a droplet 
is assumed to be uniform in all directions, then the average velocity of this 
flow in any direction is zero. Therefore, the vapor leaves the droplet 
surface with a velocity equal to that of the parent droplet (Nigmatulin, 1967; 
Solbrig and Hughes, 1975; and Jackson and Davidson, 1983). In this case no 
differences should be produced on the momentum equation of the droplets' 
Equation 2.8 from those of solid particles (Crowe, 1980). 

I 

Solbrig and Hughes (1975) tested the relative importance of the transient 
force terms in the momentum equation of solid particles under different flow 
conditions. They reached the same conclusion of many other workers (Sha and 
Soo, 1978; Soo, 1967; and Hjelmfelt and Mockros, 1966) that the Basset or the 
transient term and the virtual mass term in the momentum equations of the 
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solid particles can be neglected if they are moving in a gaseous media. 
Therefore those terras are neglected in the present study. 

The concentration equation 

To avoid the problem of density fluctuations of the carrier phase at this 
stage, only isothermal flows are considered and vaporization is assumed to be 
due to the vapor concentration gradient only. 

The concentration C is defined as the ratio of the evaporated mass within 
a control volume to the mass of the carrier phase in the same volume. The 
instantaneous, volume-averaged concentration equation for the evaporating 
material is 


k.k 


<Vi c >, t + + 1 * 4 
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where 6 is the molecular mass diffusivity of the evaporating material in air, 

k k 

The source term in Equation 2.9 ( l * m ) represents the evaporated material 

k 

from the droplets of different sizes. 


2.3 - Time-Averaged Equations 

Introduction of time-averaged quantities . For steady mean flow, the time 
averaged or mean values of U^,V^,P and are defined as the following: 
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Following common practice, all the quantities are separated into to a 


fluctuating and a time average component as follows: 


U i = U i + U i 


V i = 


V + v 
i i 


P = P +p, $ = $ + $ 


& C = C + c 
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For brevity, the overbars indicating averaged values will be dropped from all 
the quantities herein. 

The mean continuity equation of the carrier phase . Introduction of Equation 
2.11 into Equation 2.2, and subsequent time averaging yields: 
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The mean continuity equation of the k L phase . Introduction of Equation 2.11 
into Equation 2.3, and subsequent time averaging yields: 
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Tbe mean momentum equations tor the carrier phase . Introduction of Equation 
2.11 into Equation 2.5, and subsequent time averaging yields: 
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Multiply Equation 2.12 by Uj , then subtract from Equation 2.15 and rearrange, 
the result will be 
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The mean momentum equations for the k phase . Introduction of Equaiton 2.11 
into Equation 2.6, and subsequent time averaging yields: 
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Multiply Equation 2.13 by Vj , then subtract from Equation 2.17 and rearrange, 
the result will be 
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The mean concentration equations . Introduction of Equation 2.11 into Equation 
2.9, and subsequent time averaging yields: 


(Pj»jU C) 






i * k - k 

+ l $ m 
k 


- PjC^cUj + C^u^ + (^UjC 


+ 4.u.c) . 
1 J >J 


2.19 


Multiply Equation 2.12 by C, then subtract from Equation 2.19 and rearrange, 
the result will be 


ViV , 


= (.Pj «(® t C 




7 )J 


4> k m k (l-C) 


- PjC^cUj + 


C *l U j 


VjC + Vj c) ,j + p l C( *l u j ) ,j 


2.20 


In Equations 2.12 to 2.20 the overbars indicate Reynolds averaged 
correlations. 


4 > 


33 



2.4 The Problem of Closure 


In order to close the system of equations, the turbulent correlations in 
Equations 2.12 - 2.20 must be modeled in terms of the time-averaged quantities 
and some turbulence quantities that are governed by the laws prescribed by a 
"turbulence model." Examination of the literature on the mathematical models 
of two-phase flows shows that most of the existing models are based on ad hoc 
modifications of the single-phase turbulence kinetic energy and length-scale 
equations. As a result, those models fail to predict the physical behavior 
of two-phase flows. The next section, 3, describes a two-equation turbulence 
model suitable for dilute vaporizing sprays. 


9 
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3.0 A TWO-EQUATION TURBULENCE MODEL 

3.1 Introduction 

The objective of this section is to develop a general and economical 
turbulence model for free bounded dilute vaporizing sprays. 

The first task is to select the type of model that is to be employed. 
Thus, in Section 3.2, the necessity to consider a model that employs transport 
equations for both the energy and the scale of turbulence will be pointed 
out. Starting with the instantaneous two-phase momentum equations for an 
isothermal flow, the transport equations for the turbulence kinetic energy and 
its dissipation rate for the carrier phase are obtained in Section 3.3. 

Closure of the proposed set of transport equations is achieved by modeling the 
turbulent correlations up to a third order in Section 3.4. The modeled 
equations in the Cartesian tensor notations are presented in Section 3.5. 
Finally, the modeled equations in the cylindrical coordinates are presented in 
Section 3.6. 

3.2 Choice of Model Type 

As discussed in Subsection 1.2.3, the previous attempts to model the 
dilute suspensions of two-phase flow to account adequately for major exchanges 
of momentum and mass between phases has not yet been established, even for 
dilute systems containing particles smaller than the Kolmogorov length 
scale. Few investigators have tried to consider the effects of the 
particulate phase on the turbulence structure (Nagarajan and Murgatroyd, 1971; 
and Genchev and Karpuzov, 1980) but they introduced many phenomenological 
approximations and coefficients that render their schemes are applicable to 
more general flow conditions and configurations. Melville and Bray (1979) and 
Michaelides (1984) have employed the mixing length hypothesis to handle the 
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gas solid two-phase flow in free jet and fully-developed pipe flows 

respectively. This approach is limited to flows where turbulence structure 

changes at a slow rate in the main flow direction. However, the empirical 

constants involved vary from one flow situation to another and are thus valid 

for restricted flow regions only. Danon et al . (1977) employed a one equation 

model (K equation) to consider the effects of particles on the turbulence 

quantities. The deficiencies of that model to obtain accurate predictions of 

« 

two-phase turbulent jet flows necessitated that they multiply the production 
and the dissipation of K by coefficients that are dependent on particle size 
and concentration. The encouraging results obtained by Elghobashi and Abou- 
Arab (1983), Pourahmadi and Humphrey (1983), and Buckingham and Siekhaus 
(1981) suggest that higher levels of closure are required to predict shearing 
two-phase flow accurately. The model will be based on the two-equation (K-e) 
model of turbulence for single-phase flows with universal constants (Launder 
et al., 1972). 

3.3. The Exact Equations for K and e 

3.3.1 The Turbulence Kinetic Energy Equation (K) 

The equation governing the mean kinetic energy (K = 1/2 u^u^) of 
turbulence is obtained by substituting with the mean and fluctuating values 
instead of the local values in the instantaneous momentum equation of the 
carrier fluid (Equation 2.7), then multiplying by , and finally time- 
averaging. 

The resulting K equation reads: 

Convection (I) Production 
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- ff.' ♦l p ,i U i + P ,l+l"i ' ♦lV.l + "iViVm + p l*l u iVi,*) 

(II) Turbulent Diffusion 


- fU u A u. + U U A-U.) 

1 £ i^ 1 i ,£ £ i ,JT 1 1 * 

(III) Production and Transfer 


1 n , Jc .k> , k , k „ . k.k . k.k ... 
+ - l (F + m } (A U t ) + A u i (v i“ u i ) + ♦ u i( v i “ u i>) 


1 k 


(IV) Extra Dissipation (e ) 

P 


+ ( p i u i(<“i, t + u t ,i )t i ) ,i ) 4 (v i u i ( *i (u i,t + U t ,i )) ,t +V i u i ( *i (u i, lt + Vi’) .it) 


(V) Viscous Diffusion 
and Dissipation 


(VI) 


Extra Viscous Diffusion 
and Dissipation 

3.1 


3.3.2 The Turbulence Energy Dissipation Rate Equation ( e ) 

The exact equation for the dissipation rate per unit volume 

e = v,(u . u .)} is derived by differentiating the instantaneous equation 
i ^ » J 1 » J 

(2.7) with respect to x^, then multiplying throughout by v,u. and finally 

1 i >3 

time averaging. The exact equation of t thus obtained is 


f*.U e ) = f ~ 2y u . .u. (A U ) . - 2v u .u .A.U. - 2v (a u ) .u. .U. 

v 1 £ v 1 i,J i»£ 1 £ ,J 1 i,J £,J 1 i,£ 1 1 £ ,.1 i,3 i,£ 

Convection (I) Production by the Mean Motion 


- 2v.u . u (a ,U ) . - 2v,A, u u , .U. . - 2v,u .u (A.U ) . ) 

1 i , J i,£ M £ ,j 1 £ i,j i ,£ J 1 i,J £ 1 i ,£ ,J J 


% 


(II) Production by Self-Stretching of 
Vortex Tubes 
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+ [-2v,4>,u. .u. .u - 2\i <t>„u. .u. . .u ] 

1 1 i,j i,4j 4 1 V 1 i,j i,4j 4 

(III) Turbulent Diffusion 


+ [-2v.<J> 1 u. .(U U. ) . - 2v u. .U U. „ “ 2v U A u. .u. ] 

11 i,j 4 1,4 ,j 1 i,J 1,J 4 1,4 1 4 1 i,j l , 4j 

(IV) Production and Transfer 


+ — [-2v u . .(p .* ) . - 2v,u. . ( p .) . - 2v u. . ( <t> 1 P .) ,] 

P L 1 i»J ,i’l ,J 1 i.J l,i ,J 1 i,J 1 ,i ,J 

(V) Spatial Transport by Pressure (fluctuation and mean) 


-2v, 


kx k. 


+ 7 (F + m ) [u. .[4> (U.-V.)] .+u. . [(u.-v.)* ] .+u. . [ <{> (u.-v.)] .] 

i,J 1 1 ,J i,J ii ,Ji.J ii, J 

(VI) Extra Dissipation 


+ 2v r v u. : [(’i. +u . )* ] 

1 1 i,J i,4 4,1 1 , J 4 

(VII) Viscous Diffusion and Destruction 


+ [2v,\; 1 u. . [ <J> ,(U. +U .)] . +2v ,v,u. . ( 4> ,(u. +u •)) -.] 

1 1 i,j V 1 1,4 4,1 , j 4 1 1 i,j Y 1 1,4 4,1 , J 4 

(VIII) Extra Viscous Diffusion and Destruction 


3.2 


The terms in the K and e equations, (3.1) and (3.2), are classified into 

groups enclosed by large curved brackets; each group is labeled according to 

its particular contribution to the conservation of the transported quantity. 

The turbulent correlations in Equations 3.1 and 3.2 that include the 

. k . » 

fluctuation of the volume fraction, or $ , or their gradients are due to 
the presence of the particles in the same control volume with the carrier 
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\r 

phase; setting these correlations to zero, $ to unity and $ to zero' will 
reduce Equations 3.1 and 3.2 to their familiar counterpart for single-phase 
flows. 


3.4 Closure of the Proposed Set of Transport Equations 

The proposed model is restricted to high Reynolds number flows of dilute 

vaporizing sprays. Therefore, the viscous diffusion in all the governing 

equations is neglected due to its relatively small magnitude as compared with 

the turbulent diffusion. Due to the diluteness assumption of the suspension 

all fourth-order correlations containing the volume fraction fluctuations such 

as <b,u.u u, u. .u, (u *,) . and Au, .u. .u are neglected due to their 
9 1 i l M i.J rl »J i.J it£J £ 

relatively small values. The continuity equation of single-phase flow is used 
in the modeling approximations of some of the turbulent correlations. 

3.4.1 Closure of the Continuity Equation of the Carrier Phase 
The second term on the LHS of Equation 2.12 represents a mass flux 
contribution to the turbulent diffusion of the carrier phase. Following Hinze 
(1972), Melville and Bray (1979) and Elghobashi and Abou-Arab (1983), a 
gradient-type diffusion is assumed for this correlation given by: 


+ i u i ‘ ' E f 4 i,i 


3.3 


where e ^ is the mass eddy diffusivity of the turbulent flow of the carrier 
phase. This quantity can be related to the momentum eddy diffusivity (v t ) as 
the following: 



3.4 


39 



where o c is the effective Schmidt number for the carrier phase. It may be 
expected to be constant of value 0.7 in line with the average levels of 
effective Schmidt number reported for a number of free shear flows (Launder, 
1976; and Spalding, 1971). 

The momentum eddy diffusivity of the carrier phase is related to fluid 
kinetic energy (K) and the rate of dissipation (e) of K by: 

v = c K^/e 3.5 

t p 

The value of c in general is 0.09 but it can be a function of suitable flow 
P 

parameters to extend the range of applicability of the K-e model. For 
example, in axisymmetric jet flows which are considered in this work, those 
parameters are the deceleration of the velocity at the axis of the jet 
(Equation 3.36) and the jet width (Equation 3.37). 

Corrsin (1974) discussed the limitations of the simple gradient 
hypothesis for modeling turbulent diffusion in turbulence. He pointed out 
that this model may lead to inexact results if the size of the 
energy-containing eddies is much smaller than the distance over which the 
gradient of the considered quantity varies appreciably. 

Lumley (1975) tried to overcome this problem by proposing a model for the 
turbulent flux of passive scalar in inhomogenous flows. But since Lumley' s 
model is not well tested yet, the simple gradient hypothesis will be used in 
the present work due to its fruitful results in many types of flows (Lunder et 
al., 1972). 

t h 

3.4.2 Closure of the Continuity Equation of the k Phase 

• Ic Ic • 

Similar to Equation 3.3 the correlation <j> v^ on the LHS of Equation 2.13 
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is modeled as 


k k k k 

♦ v i“ ’ e h‘,i 


3.6 


is the mass eddy diffusivity of the turbulent flow of the droplets. An 

expression to obtain this quantity in terms of its counterpart for the carrier 

phase and the droplets’ Schmidt number a is developed in section 5. The 

P 

droplets' Schmidt number is defined as 


k k , 

% = G h /e f 


3.7 


3.4.3 Closure of the Momentum Equations of the Carrier Phase 
Here, the modeling of the turbulent correlations needed to close the 
momentum equations (2.16) are presented. The correlations of two scalars 
containing the volume fraction in the momentum and concentration equations 


such as (j> c, 4 , p or $ cu. are neglected. This approximation is based on the 
l l ^ J 

following: 1 ) the lack of understanding of the nature of those correlations, 
thus the modeling which is supported by the experimental data (Lumley, 1978b) 
is not available, and 2 ) their relatively small values compared to the 
turbulent diffusion terms (Buckingham and Siekhaus, 1981; and Launder, 1976). 


The correlation u^u. in Equation 2.16 represents the transfer of momentum 
by the turbulent motion. The oldest proposal for modeling this correlation is 
that of Boussinesq (cited by Launder, 1976): 


U i U j = 


v (U. .+ U. .) + | 6, .K 
t 1,3 3,1 3 i 3 


3.8 


The term involving the Kronecker delta, 6 , is necessary to make the 
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expression applicable also to normal stresses (when i = j). The expression 
(3.8) has been severely criticized by some workers, and it should not be used 
without caution. The use of that expression is justified on the basis of an 
approximate local equilibrium. If the addition of the droplets causes the 
turbulence of the carrier phase to adjust more slowly to the mean velocity 
field, or if it introduced additional mechanisms for generation of turbulent 
energy, the expression (3.8) will be a poor approximation. Melville and Bray 
(1979) argued that neither of these effects will happen if the mass loading 
ratio is less than unity. The other approach is to solve a transport equation 
for uTuT, which in turn contains. higher order correlations that require 
modeling. To be consistent with the present level of closure the present 
study will use (Equation 3.8). 


The correlation (j> u.u. will be modeled by adapting Launder's proposal 
i r 3 

(1976) that gives 


4>.u.u. = - c (K/e) [u.u(u .<►, ) . + u .u (u . <J», ) 1 

1 l j x l j 1 , !■ j l l 1 ,i 


3.9 


where c is a constant of value 0.1. 

• k Ic • 

The correlation $ (u^-v^), which appears in the momentum equations of 

both the carrier phase and the droplets, can be written as 


.k/ k-, _ k k k 

b vu.-v.) = * u. - * v. 

XL X X 


3.10 


The second term on the RHS of Equation 3.10 was discussed previously. The 
first term is modeled as 


k 

0 u 




3.11 
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The last correlation to be modeled in the momentum equations of the carrier 


phase is .u. . , which appears multiplied by . , and will therefore be 
1 i > J 1 

neglected due to its relatively small value. 

3.4.4 Closure of the Momentum Equations of the k c ^ Phase 


Similar to the carrier phase treatment the correlations <$> p . and 



k k k 

(v + v ) in Equation 2.18 are neglected. The two correlations 
i » J J > i 

k k k k 

<P v^ and i p (u^ - v^) have already been discussed. The only two correlations 

k k k k k 

still to be modeled are those of the forms v.v. and v.v.. Similar to the 

i J i ] 

k k 

carrier phase, the correlation v ^ v j is modeled as 


k k 

V i V j = 


- v k (V k . +V k ) + \ 6 .K k 
P i,J J,i 3 i] p 


3.12 


where 


„k 1 k k 
K p ' 2 v l v i 


3.13 


The momentum eddy diffusivity of the droplets in the size range k ( v p ) is 
related to its counterpart for the mass as the following: 




3.14 


where is a coefficient of value 0.7 as given by Melville and Bray (1979) 

k k k k 

Again, is determined in section 5. The correlation $ v ^ v j is modeled 

similar to Equation 3.9 as 


,k k k /•„/„%( k k, k,k. , k k, k.k. , 
* v i v j = c ^ K/£ ) K v * ( V >,A + v j v A (v i^ 


3.15 


43 



3.4.5 Closure of the Vapor Concentration Equation 


As discussed previously the correlations 4>^c j , ^c and in Equation 

2.20 are neglected compared with cuT or <j>^Uj . An investigation of the 


behavior of cu^ suggests that it can be modeled similar to Equation 3.3 as 


cu. = - e,C . 
J f »j 


3.16 


3.4.6 Closure of the Turbulence Kinetic Energy Equation 
The exact equation of the turbulence kinetic energy K for the carrier 
fluid is given by Equation 3.1. The terms are grouped according to their 
physical contribution to the conservation of K. The modeling of the turbulent 
correlations appearing in the K equation (14 correlations) are presented in 


this section. (J^u^u^u^ ^ is neglected since it is a fourth order 

- - - . ■ - 

correlation. The correlations u.u., <i>, u. , d>,u.u„and <t> u. in Equation 3.1 were 

i Z’ 1 i T 1 i * i 

discussed in subsection 3.4.2. 


The pressure diffusion terms (u^p ^ and ^p in Equation 3.1 are 

neglected, following the Imperial College group and the recommendations of 
Hinze (1975), and because very little is known about it (Launder, 1976). 


The turbulent diffusion correlation u^u^u^ ^ can be written as 


[u a ( 1/2 uj)]^ 


3.17 


which modeled as 


V°k )K ,* ] ,Z 


3.18 


where is an empirical diffusion constant of order one. 
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The correlation ^ can be written as [$j(l/2u^] ^ , which can be 


neglected due to its relatively small value compared with [u (l/2uf0] .• 

X 1 ) iy 


Following Elghobashi and Abou-Arab (1983), the correlation u^(v*-u^) is 


modeled as 


u i^ v i“ u i ^ = “ K U - J* ((^j- dt0 J 


3.19 


where ou is the harmonic frequency of turbulence and E(u>) is the Lagrangian 
energy spectrum function of the carrier phase. anc * S1 R are functions 

of the carrier and dispersed phases properties, the droplet diameter, and the 
harmonic frequency. They are discussed in detail in subsection 4.1. 


Ic k 

The correlation $ u^v^-u^) is modeled as 


'f'Vj^v^- u i ) 


- u i )u i l<u ** k \» + (<v r u * Hk) ,*J 




3.20 


The extra viscous diffusion and dissipation, two terms in group VI in 
Equation 3.1, are neglected due to their relatively small magnitudes as 
compared with the other similar terms (see Daly and Harlow, 1981, and Launder 
et al., 1976; Launder et al. 1975). 

Neglecting the viscous diffusion, the last correlation to be modeled in 
the K equation, UjU^^u^^ + ^ represents the dissipation rate of 

K. When the local Reynolds number of turbulence is large the dissipative 
motions can be assumed to be isotropic, therefore 


S^ (u i ,* + V 


3.21 
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3.4.7 Closure of the Turbulence Energy Dissipation Rate Equation 

The exact equation of the dissipation rate of turbulence energy, e, for 
the carrier fluid is given by Equation 3.2. The terms are also grouped 
similar to the K equation. 

Tennekes and Lumley (1972) have shown, based on an order-of -magnitude 
analysis, that the terms involving mean strain rates in Equation 3.2, groups 
(I) and (IV), are negligible at high turbulent Reynolds number compared with 
the production by self-stretching of vortex tubes. Therefore, groups (I) and 
(IV) are neglected in this study. 

The correlation 2v^4jU^ jU^ £j u £» which accounts for the diffusion 
of e by velocity fluctuations, is handled as the following: 


2 Vi,j U i,*j U A = ( Vi,j U i,j U * ) ,A 

" (e u *\* 

- - i(v t /o e )e,tJ, t ' 3.22 

Group (V) represents the diffusional transport of e by pressure fluctuations. 
Also, it contains a term that represents a transfer due to the mean pressure 
gradient, 2vu. ,(<J>,P .) .. Following Rodi (1971), the present study neglects 
the pressure diffusion terms (group V). 

The viscous diffusion terms, part of group (VII) and group (VIII), are 
neglected due to their relatively small values compared with the turbulent 
diffusion. 

The first term in group (II) expresses the generation rate of vorticity 
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fluctuations through the self-stretching action of turbulence. Rodi (1971), 
and Hanjalic and Launder (1972) have argued that this term should be 
considered in conjunction with group (VII), representing the decay of the 
dissipation rate ultimately through the action of viscosity. At a high 
turbulent Reynolds number, these two terms are of opposite sign, however, and 
their difference necessarily remains finite. Following Rodi (1971), the terms 
are collectively approximated as follows: 


- 2 v ,u 




(c el I V e ~ c e2*P K 


3.23 


where is the total production of K [group (I) and (III) in Equation 3.1], 
and c £l and c are constants of value 1.43 and 1.92 respectively. 

The last correlations to be modeled in e equation are those of group 
(VI), which represent the extra dissipation due to the relative velocity 
between the phases. They are modeled collectively as one term which is given 
by 

- C £3 E p £/K 


where (term IV in Equation 3.1) is the extra dissipation of K. The 
constant c £ ^ was optimized by Elghobashi et al . (1983) for a two-phase jet 
flow. The value of this constant is 1.2. 


3.5 Modeled Transport Equations in the Cartesian Tensor Notations 
Using the modeling approximations discussed in subsection 3.4, the 
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transport equations in the modeled form are obtained. Those equations are 
given in Appendix B. 


3.6 Modeled Transport Equations in the Cylindrical Coordinates 
The flows considered in this work are 
a) axisymmetric (without swirl), 

1 A 

Vg = Ug = 0, gg (time-averaged quantites) = 0, 3.24 


b) of the boundary-layer type. 


V k > V k , U > U 
z r z r 



3.25 


The present study adopts the notations commonly used for the boundary layer 

flows: z, r, 0 for the coordinates; U z , V r , V Q for the velocity components of 

k k k 

the carrier phase; , V^, Vg for the velocity components of the droplets of 
class k. 

Using expressions 3.24 and 3.25, the modeled transport equations for the 
mean and turbulent quantities presented in Appendix B can be expressed in 
cylindrical polar form. This can be done in a straightforward manner by the 
methods of tensor calculus as exposed in Synge and Schild (1978), for 
instance. 

The mean continutiy equation of the carrier phase is 




+ — (r$,U ) 
r 1 r ,r 


1 o 1 , z 
c 


P 1 v t 

- -- (r-- ® ) 

z r a 1 ,r . 
c 


= l m k $ k 3.26 

k 
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The mean continuity equation of the k tri group is 


„ k . V 2 . ,k k\ , p .k . v 2 , p .k . .k_k 

Po(* v J + (r* V ) - P 0 (--- * J - -- ( r--- $ ) = - m $ 


z' ,z r 


r ,r 2a ,z ,z r 
c 


,r' ,r 


The mean global continuity is 


* , + l * = 1 

1 k 


The mean momentum equation for the carrier phase in the axial (z) 


direction 


P^U U + p_4>U U 
1 1 z z,z 1 1 r z,r 


$,P -l 4 k (F k -H?i k )(U -V*) 
1 »z it z z 


i V t i k v t 

+ -(♦, P.rv U ) + p.u (-- *, ) + c.p, - I- rvu ) 1 

r 1 1 t z,r ,r 1 z,r v o l,r' <j> 1 r '■e t z,r l,r ,r',i 


The mean momentum equation for the carrier phase in the radial (r) 
direction is 


P.fc.u U + p.*,U U = - 4 P 
1 1 z r,z K 1 1 r r , r 1 ,r 


l $ k (F k +m k )(U -V k ) 
k r r 


- 1 ^ - l * jE Cl-<$>* r 

k c ’ k c 


V t P 1 V t 2 P 1 

+ P -- *, U + ~ L (-- rU *, ) - | ~ (rK*,) 

1 o l,r r,r r a r l,r ,r 3 r 1 ,r 

c c 
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3.30 


P 1 ,4 K 2 \ . . 

+ ST ( 3 r r <T \r\r\r 


The mean momentum equation tor the k ck phase in the axial (z) directon is 


■+ = - * k p + F k * k (U - v k ) 

2 z z,z 2 r z,r ,z _ 


z z 


+ -( $ k r p o W ) + p.V k (v k $ k ) + c p — (-r v k V k ( v k $ k ) J 

r 2 v p z,r ,r 2 z,r P ,r $ 2 r '■e p z,r P ,r ,r',r 


+ ( p 2 - 


3.31 


The mean momentum equation for the k tk phase in the radial (r) 
direction is 


.k „k k ,k k k ,k_ - „k k.„ tI k. 

P-* v v r , + P 0 $ V V --*P+F* (U -V ) 
2 zr,z 2 rr,r ,r rr 




+ m k -E 4> k + F k -- (l-^)* 1 


°c ’ r 


, k 

P' .r 


„ -P- xk_„k P 2 , p ..k.k v. 2 P 2 , „k k> 

2 a ,r r,r r v o r ,r\r 3 r p ,r 
c c 


O V 

2,4 K k, t 


+ c . — (i r - K (-- ♦. . ) 1 

<(> r ^3 e p v o 1 ,r ,r ' ,r 


3.32 


The mean concentration equation is 


p,$,U C + p,$, U C = - (p.r*. -- C ) 

1 1 z ,z 1 1 r ,r r 1 1 ,r,r 


+ z • A (1-C) + Pl c tr (^ V>* 


3.33 
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The turbulence kinetic energy equation (K) is 


4 V t V t 

°lVz K .z * p lVr\r ■ 'iVA.A,, - 3 ’filial *l,r>,r D r,r 
+ «Wl>VV Vr>,r u z,r"z,r ' * < " )d “ ) 

C K 4 


+ l (F k t Jl k )[(U r -V k )(v p k » i “) - c t( :i)(v p k * k r ) >r (l- j o “ (i-fS)E(«)4»)] 


+ ? <p l*l rK ,r\r - p lV > 


3.34 


The turbulence energy dissipation rate equation ( e) is 


c . v v 


P,*,U e + p * u e = c , l p,*, [ v U U - \ ) U 

llz,z 1 1 r ,r el K lit z,r z,r 3 4^ c^ o c l,r ,r r,r 


ft,- 


+ r <1> \ ♦l,r>,r D .,r l, .,rl ‘ p z3 I (fV ^ 

1 C K 2. 


- (U - V k )(v k * k ) + c .(--)( v k $ k ) (1- J )E(u)dw)] 

r r P ,r' <p 'c 7 p > r > r 0 


+ t (p l r *l ®; e ,r> ,r - c « p lV £ K • 


3.35 


The constants of the turbulence model 


The constants in Equations 3.26 to 3.35 are two sets: one is identical to 
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that of the single-phase K-e model and the other belongs to the two-phase 
model. The former is well established and it was not changed here. The 
latter is new (c^, o,) and has been obtained from one set of data (Elghobashi 
et al . , 1984). The coefficients of the single-phase model and the optimized 
values of the new coefficients are given in Table 3.1. 

Table 3.1 Coefficients of the Turbulence Model 


o 

c 

°k c u °e c * 

c el c e2 c e3 °v 

0.7 

1 . Eq.( 3.36) 1.3 0.1 

1.44 Eq . ( 3. 38) 1.2 0.7 


c y = 0.09 - 0.04f ; 


3.36 



dU 


z,c 


0.5R 


dz 


z ,c 


dU. 


z,c 


0.2 


dz 


U 


z 


3.37 


n „ and U „ are the axial velocities of the fluid at the jet centerline and 
z , c z , 00 

the ambient stream respectively, R is the local jet width (Launder et al.., 
1972). 

The constant c e2 Table 3.1 is given by 

c e2 = 1.92 - 0.0067 f 3.38 

The quantities ^ > 

terms in Equations 3.34 and 3.35 will be discussed in subsection 4.1. The 

• k 

evaporated mass per unit time and unit droplet volume (m ) is calculated in 


and E( u) , which are involved in the integration 
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A. 3. The interphase friction coefficient (F^) is evaluated in subsection 
4.4. The last quantity in Equations 3.26 to 3.35 that should be calculated is 
the momentum eddy diffusivity of the droplets (v ). This quantity is 
calculated in section 5. 
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4.0 SINGLE PARTICLE BEHAVIOR IN A TURBULENT FLOW 

In this section some turbulent correlations that are needed in the 
turbulence model closure are obtained. First the equation of motion of a 
single particle in a turbulent flow will be discussed, then Chao's (1964) 
solution for that equation will be presented, and finally the turbulent 
correlation between the fluid velocity and the relative slip velocity will be 
obtained. 

4.1 Transport Behavior of a Single Particle 

Tchen (1947) extended the Basset-Boussinesq-Oseen (BBO) equation for the 
unsteady Stokes motion of a spherical particle in a stagnant fluid to that of 
a particle in a moving fluid that reads: 


^3 . 4 TTa . 2 TT3 f • . v / f v 

3 — P 2 V = "3 — p l u 3 P p v ~ u ' “ 6*^3 (v - u) 


4.1 

In Equation 4.1 u(t) is the velocity of the fluid in the neighborhood of the 
particle, but far enough to be unaffected by it; v(t) is the velocity of the 
particle; p^and p ^ are, respectively, the density of fluid and the density of 
the particle; a is the particle radius; g is the acceleration of gravity; and 
p^and are, respectively, the dynamic and kinematic viscosities of the 
fluid. The dot denotes ( 4 — + v 4 —). The physical meaning of the different 
terms of Equation 4.1 are discussed by many workers (Fuchs, 1964; Hinze, 
1975). As a summary, the first term on the right hand side (RHS) of Equation 


a ,t *<'l> - «‘l> 
- 6irp,a . l m 

/ivj /t"-"t“ 


4ifa 


dt. 


g (P 2 - P x ) 
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4.1 is added to BBO equation on the basis of intuitive consideration by Tchen 
(1947) to account for the unsteady motion of the fluid or its pressure 
gradient. The second term is the inertia force due to relative acceleration 
of the virtual mass attached to the particle; the third term is the Stokesian 
drag; the fourth term is often referred to as the "Basset" force, which 
results from the relative acceleration between the particle and the fluid; and 
the last term is the gravity force. The importance of each term in Equation 
4.1 under different flow conditions and the various approximations for 
handling that equation is explored by Hjelmfelt and Mockros (1966). They 
found that for high density ratios (e.g., those greater than 10 3 ) all the 
terms on the RHS of Equation 4.1 are unimportant except for the drag term. 

For the sake of generality, however, all the terms will be retained in the 

r 

present analysis. 

Corrsin and Lumley (1956) argued that the first term on the RHS side of 
Equation 4.1, which Tchen included to represent the force created by the 
pressure gradient in the fluid, should be evaluated via the full Navier-Stokes 
equation. Accordingly they proposed that Equation 4.1 should be replaced by: 


4*a 3 dv i 
3 P 2 dt 
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g (P 2 - pj) 


4.2 


4ira' 


Maxey and Riley (1983) derived an equation of motion for a solid particle In a 
turbulent flow. They considered the Faxen terms to account for the unsteady 
effects on Stokes drag law. For particles of size smaller than the 
Kolmogoroff length scale, b, Hinze (1975) argued that the following two 
conditions can be satisfied: 


a) 


2 3u 


i 

~5x7 


« l 


b) 


v 1 3u 1 /3x i 
V 1 S^/Sx 


2 

i 


» 1 


4.3 


4.4 


Thus, if the conditions 4.3 and 4.4 are satisfied, then the viscous stress 
term -in Equation 4.2 and the Faxen terms on Maxey and Riley's equation can be 
neglected. Therefore, in the case of particles with a diameter less than the 
Kolmogoroff length scale, all the equations of Corrsin and Lumley, and Maxey 
and Riley, become identical to Equation 4.1. 

Chao (1964) considered Equation 4.1 with the two„res trict ions 4.3 and 4.4. 
to obtain the connection between some transport properties of a particle and 
those of its surrounding fluid. Neglecting the gravity force and applying 
such restrictions. Equation 4.1 reads as follows: 


„ , , , (v. - u! )dt , 

- M i + - ”i> + -Vrp- A ’ 0 


4.5 


in which a = 3'^/a‘ i 


4.6 


3P, 


8 = 


2P 2 + P 1 


4.7 
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Chao (1964) applied the Fourier transform of the velocity component u^(t) 
as defined by 

( “) = J l^Ct) exp (-imt)dt 4.8 

and similarly for v^. Unlike u^, is not only a function of oi, but also a 
function of the physical parameters a and 3. 

The Fourier transform solution of Equation 4.8 with t Q = - <*> gives 


_ + 1 I«+^] _ 


V 


, , 3 am , , r oj . 3 aai . 

[a + ✓ -g- ] +![- + ✓ — ] 


u. 

l 


4.9 


By introducing the energy spectrum function E(w), a relation is obtained 
between the intensity of the particle turbulent motion and that of the 
surrounding fluid: 


v . ^ 

— = I -J- E(«)dw 
— 2 o « 2 

u i 


where ^ - (£) 2 + VS (%) JU + 3(J + ^6 (^) 1/2 + 1 ; 


V e 2 0 2 + ^ e_1 0 3/2 + 3 0 + ' 6 0 1/2 + 1 


4.10 


4.11 


4.12 


The Lagrangian frequency function E(o>) is in general affected by the 
presence of the particles. In the low frequency range (inertial subrange), 
the modulation of the Lagrangian frequency function of the carrier fluid by 
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the particles can be neglected (Al-Taweel and Landau, 1977). Thus, in the 
present work the Lagrangian frequency function is given by (Hinze, 1975) 

E(o>) = (|) ( , 4.13 

1 + 0) T L 
Ll 

where <o ranges from 1 to 10^ (sec - ^). 

The local Lagrangian integral time scale t is calculated from (Mostafa 

Li 

and Elghobashi 1984b): 

t l = 0.233 K/e 4.14 

Chao proceeded his solution by defining a relative velocity u)^ between 
the particle and the local gas velocity as 



4.15 


This value, when substituted into Equation 4.1, again with t = - °°, followed 
by the Fourier transformation and solving for , gives 


<v 



- (1 - 8)iu> 

[o8 + $(~) V ‘] + i[u + 



4.16 


The first term in the denominater of Equation 4.16, a0, was written wrongly 
(as a only) by Chao (1964) and Soo (1967). Hetsroni and Sokolov (1971) 
handled the incorrect form of Equation 4.16 to study the effect of the 
dispersed phase on the fluid turbulence energy spectrum. Using that erroneous 
equation Hetsroni and Sokolov obtained Equation 29 (in their article), which 
relates the turbulence intensity of a two-phase flow to that of a single 
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phase. If Equation 4.16 is used in the analysis of Hetsroni and Sokolov, 
their Equation 29 will give equal turbulence intensities for the single- and 
two-phase flows. Although Hetsroni and Sokolov's theoretical analysis gives 
good qualitative results, it has no basis and therefore should not be taken 
seriously. 

Chao (1964) obtained the ratio between mean square relative velocity to 
that of the surrounding fluid: 

u 2 ^ fl 

— = l“ ** E(io)da> 4.17 

-r ° *2 

u . 

1 

where = [ (1-8) ui/ a0] 2 4.18 

This solution applies to dilute suspensions, where there is no interactions 
between the particles. 

Elghobashi and Abou-Arab (1983) considered Chao's solution to get the 

turbulent correlation u.(u. -v.) that is needed in their turbulence model. 

i i i ' 

This quantity is given by 

_____ fl - o 

u i( u i " v i^ = “ | u i ^ ~ Jo — H E(w)dw) 4.19 

This correlation will be used to close the time-mean equations in the present 
work. 

4.2 Droplet Shape 

4.2.1 Theoretical Analysis 

It is well known that the change in droplet shape affects not only the 
interfacial area and drag force but also the evaporation rate. Most of the 
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theoretical studies (Krestein, 1983; O'Rourke, 1981; Mostafa and Elghobashi, 
1983) for liquid droplet-gas flow have been based on the assumption of 
spherical droplets. This assumption must be justified, especially if the 
droplet is suspended in a turbulent flow. 

The shape of a liquid droplet moving in a continuous phase is determined 
by the forces acting along the surface of the droplet. At any time the net 
force is the balance of the pressure, gravity, buoyancy, drag, and inertial 
forces of the exterior fluid. At the fluid-fluid interface there will be an 
equilibrium of normal forces. The forces acting inward are due to the dynamic 
stresses and static head of the exterior fluid and interfacial tension. Those 
acting outward are due to the dynamic stress and static head of the interior 
fluid. If the droplet is spherical, all the forces will lie on a radial line 
and the . interfacial tension force will be the same on all parts of the 
surface. 

For a liquid droplet moving in a gaseous flow, a study is presented here 
of the physical factors that might be expected to control the spherical 
shape. Those factors are as follows: 

1. Surface "tension" -- This force is a consequence of the net inward 

attraction exerted on the surface molecules by those which are lying 
deeper within the droplet over the prevailing force in the gas 
outside. This increment in total pressure, across the interface, Ap, at 
a certain point on the droplet surface is given, in general, by 

AP = Y(l/a 1 + l/a 2 ) 4.20 

where Y is the surface tension of the liquid-air interface (N/m). 

In the special case of a spherical droplet, a i = a 2 = a > and then 
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AP = 2y/a 


4.21 


2. Internal hydrostatic pressure — There exists within the droplet a 
vertical pressure gradient of exactly the sort found in any mass of 
fluid at rest in a gravitational field. In the limit of large droplets, 
the difference in hydrostatic pressure between top and bottom of a 
droplet (2p^ga) becomes quite important in controlling droplet shape. 

3. The relative velocity between the droplet and the gaseous phase — The 

fluid dynamic pressure exerted because of the relative velocity between 

the gaseous phase and the droplet tends to cause a distortion in the 

2 

spherical shape. The effect of this inertial force (l/2p^(U-V) ) 
increases as the Reynolds number does. As the Reynolds number is 
increased, droplet oscillation (unsteady state distortion in shape) will 
set in; ultimately as Re increases droplet breakup will occur. 

4. Internal circulation — Due to the vortical motion, there is a 
centrifugal reaction that varies as the square of the circulation 
velocity (Oliver and Chung, 1982). Many workers (Beard, 1976; 

Pruppacher and Beard, 1970) reported that in the case of a water droplet 
(of diameters less than 1 mm) falling at terminal velocity, internal 
circulation has negligible effect on distortion. 

After stating the factors that might be expected to control the shape of 
a liquid droplet suspended in a moving gas, a study of the order of magnitude 
of each factor is necessary. The equation for the shape was given by 
Laplace's formula (McDonald, 1954) for the mechanical equilibrium of an 
interfacial surface, which can be written at the droplet's equator as: 

s* 
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4.22 


+ ■ H ♦ 'Vd** + 2 1(> 1 (,, - V > 2 


The first terra on the RHS of Equation 4.22 is the spherical curvature 
pressure, the second is the hydraulic head, and the third is the differential 
dynamic pressure. From the above equation it can be seen that as long as the 
spherical pressure (2y/a) is dominant, a variation in surface tension is 
unimportant. In Table 4.1 the comparative values of the three terms in 
Equation 4.22 are calculated for 100 y methanol and Freon-11 droplets (the 
materials used in the present study) moving in air with a relative velocity of 
5 m/s. Also, the dimensionless groups; defined by Equations 4.23a, 4.24 and 
4.26 are given in Table 4.1. The physical corresponding properties of 
methanol and Freon-11 are given in Appendix A. 

Table 4.1 The Forces Acting Along the Surface of 
Liquid Droplets and the Corresponding 
Dimensionless Groups. 






Re = 

We = 

Et = 


2y 

a 

N/m^ 

( P 2 ~P 1 )ga 
N/m^ 

Ip^u-v) 2 

N/m^ 

Pj (U-V)d 
y i 

Pl (u-v) 2 d 

Y 

g ( P 2 -P t )d 2 

Y 

Methanol 

872 

0.4 

14.73 

32.72 

0.135 

0.0036 

Freon- 1 1 

300 

0.744 

14.73 

32.72 

0.39 

0.02 


Table 4.1 shows that, for liquid droplets of diameters less than 
100 y and for the two different materials, the surface-pressure increment is 
large compared with the minute hydrostatic pressure differences within the 
droplet, or compared with the small aerodynamic pressures. Hence, these 
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liquid droplets do simply assume the shape implying minimum surface free- 
energy, thus accounting for their well-known spherical form. 

Because of the need to answer the question of whether internal 
circulations should be important in the droplet shape problem, the upper limit 
to the centripetal force acting on the droplet surface due to the internal 
circulation will be estimated here. Pruppacher and Beard (1970) reported that 
for a water droplet falling at terminal velocity, the velocity at the equator 
was found to be about 1/100 of the droplet’s terminal velocity. For a 
relative velocity of 5 m/s, the maximum velocity at the equator is about 0.05 
m/s. The centripetal force per unit area acting on the small lamina of the 
internal boundary layer can be calculated. This lamina of radial distance can 
be assumed to be equal to one tenth of the distance from the droplet surface 
to the internal stagnation point. Thus, for a droplet of 100 pm diameter, 
this distance is about 1.5 urn. Since the centrifugal reaction resulting from 
the vortical motion varies as the square of the circulation velocity, the 
centripetal force per unit area acting on the internal lamina of the droplet 
is found to be about 0.8 N/m^. Thus, the order of magnitude of the 
centrifugal pressure gradient appears to warrant neglect of internal 
circulation compared to the other parameters. 

4.2.2 Experimental Observations 

Many researchers have investigated the different factors affecting the 
droplet shape and the flow field inside and outside liquid droplets moving in 
a continuous phase (gas or other liquid). Wellek et al . (1966) investigated 

the effects of various properties, droplet size, and droplet velocity on 

% 

droplet shape for forty-five dispersed liquid drbplets falling through 
stationary liquid continuous phases. The maximum ratio between the viscosity 
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of the dispersed and continuous phases was that of ethylene glycol (droplet) 
and hexane (liquid system). This ratio was about 47. Empirical relations 
involving the Weber number. We, Eotvos number, Et, and viscosity ratio were 
obtained. These relations enable the prediction of the eccentricity of 
nonoscillating droplets over a wide range of Reynolds numbers (6.0 to 
1354). This number is given by 

p. (U-V)d 

Re = -J 4.23a 

M 1 

One of their relations is: 


R = -- = 1.0 + 0.091 We 0 * 95 4.23b 

a 2 


where the Weber number is given by 


P, (U-V) 2 d 

We = — 4.24 

Y 

a l _ 

and R = is the ratio of the length of the minor to the major axis of the 

3 2 

droplet. From Table 4.1 for a 100 U Freon-11 droplet moving with a relative 
velocity 5 m/s in air at the atmospheric conditions (Re = 32.73 and We = 
0.39), R = 1.037, which could be assumed a spherical shape. Garner and Lane 
(1959) measured the distortion of liquid droplets falling in gases. They 
reported the following linear relationship: 

R = 1 + 0.13 Et 4.25 


where the Eotvos E is given by 
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Et 


4.26 


g ( p 2 ~P 1 )<1 2 
~ 

For the Freon-11 droplet with conditions summarized in Table 4.1, R is equal 
to 1.003. Winnikow and Chao (1965) investigated the behavior of droplets 
falling in water at Reynolds numbers ranging from 138-971. Their conclusion 
about the deformation of nonoscillating droplets was that the droplets are 
spherical up to Et = 0.2 Since Equation 4.23 gives the distortion ratio in 
terms of the Weber number which is a function of the aerodynamic pressure, 
Equation 4.23 is recommended for the calculation of the spherical shape limits 
for present droplets in a moving gas. Pruppacher and Beard (1970) studied the 
deformation of water droplets falling at terminal velocity in air of 20°C at 
sea level pressure, and nearly water saturated by a wind tunnel means. They 
concluded that droplets with an equivalent diameter smaller than 280 pm 
equivalent to Reynolds number Re = 25 had no detectable deformation from 
spherical shape. Droplets of sizes d < 400 pm Re < 200 were slightly deformed 
into an oblate spheroid (R = 1.02). 

Beard (1976) also analyzed all the available theoretical and experimental 
data on droplets falling in gases to derive a reliable method for obtaining 
the terminal velocity and shape of a water droplet at any level in the 
atmosphere. He concluded that droplets with diameters < 1 mm (Re < 300) are 
essentially spherical. He also reported that the effects of a varying surface 

tension and internal viscosity were shown to have a negligible influence on 

\ 

the shape and terminal velocity of a falling droplet of diameters up to 1 mm. 

Now, it is clear that the assumption of a spherical shape for methanol or 
Freon-11 droplets of diameters up to 100 pm moving with a relative velocity up 
to 5 m/s (or Re = 32.72) is realistic based on the previous force analysis and 
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the experimental evidence of the other researchers. 

4.3 Mass Transfer 

There are many aerodynamic parameters that dominate the process of 
evaporation of a spray injected into a moving airstream. Relative velocity 
and free stream turbulence are the most important parameters. 

First, the evaporation of a spherical droplet, motionless relative to an 
infinite, uniform medium is considered in section 4.3.1. Then the evaporation 
rate of a moving droplet in a gaseous medium is covered in section 4.3.2. 
Finally, the effect of free stream turbulence on the evaporation rate is 
discussed in subsection 4.5. 

4.3.1 Quasi-Stationary Evaporation of Droplets Motionless Relative 
to Media 

One of the earliest investigations of evaporation in stagnant gases was 
made in 1877 by Maxwell (cited by Fuchs, 1959) who solved the steady-state 
conservation equations of mass and energy in the gas phase under the following 
assumptions: 

1. Spherical droplet 

2. Incompressible droplet fluid and surrounding fluid 

3. Spherical symmetry: forced and natural convection are neglected. This 
reduces the analysis to one dimension. 

4. The droplet fluid is of a single component. 

5. Constant pressure process 

6. Both droplet fluid and surrounding fluid are mutually immiscible, and 
there is no chemical reaction. 
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7. No spray effects: the droplet is Isolated and Immersed in an infinite 
environment . 

8. The system involves only purified fluids (there is no surface-active 
material ) . 

9. Diffusion being rate-controlling: the rate of evaporation is completely 
determined by the rate of diffusion of the vapor in the medium. 

10. Constant and uniform droplet temperature: this implies that there is 
no droplet cooling or heating. 

11. Constant gas phase transport properties. 

12. Saturation vapor pressure at droplet surface: this is based on the 
assumption that the phase-change process between liquid and vapor 
occurs at a rate much faster than those for gas-phase transport. 
Therefore, the vapor at the surface is produced at the saturation 
pressure corresponding to the droplet surface temperature Tg. 

In the case of stationary evaporation, the rate of diffusion of the vapor 
of the droplet across any spherical surface with radius r and concentric with 
the droplet is constant and expressed by the equation: 

2 dC r 

I Q = 4Rr P x 5 _ Kg/s 4.27 

where C r is the concentration of the evaporated material at radius r and 6 is 
the diffusivity of that material. Integrate Equation 4.27 with the following 
boundary conditions: 

C = C at r = 00 4.28 

r 

C = C T at r = a 4.29 

r L 
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This gives Maxwell's equation: 


I = 4IIa 6 (C - C T ) 4.30 

o L 

Strictly speaking, the evaporation of a droplet cannot be a stationary process 
since the radius and hence the rate of evaporation is constantly decreasing. 
Fuchs (1959) pointed out that when C^Pj << P2> the evaporation can be regarded 
as quasi-stationary; i.e., one can assume that the rate of evaporation at a 
given moment is expressed by Equation 4.30. Since 


I 

o 


dm 

dtT 


4.31 


3 

where t is the time and m = 4/3lIa is the mass of the droplet, Equation 4.30 
can be rewritten in the form: 


or 




dS 

dt 



4.32 


4.33 


2 

where S = 4Ha is the surface area of the droplet. Integration of these 
equations gives: 



4.34 


S 

o 


s = in§. (c - c )t 

Po L 


4.35 
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where a Q and S Q are the initial radius and surface of the droplet. The 
surface of the droplet is therefore a linear function of time. 

The assumption of equilibrium between the liquid and its vapor at the 
droplet interface suggests that the diffusive resistance of the gas to 
evaporation is large compared with the inter facial resistance. This is 
considered to be a good assumption under all conditions, except at very low 
gas pressures or for droplets whose size is in the order of the mean free path 
of the gas molecules (about 0.1 micron). Thus, equilibrium has been assumed 
in most analysis dealing with evaporation of droplets. 

4.3.2 Influence of the Stefan Flow on the Rate of Evaporation 

Observe that Equation 4.30 was derived neglecting the radial convection 
transport due to bulk flow of the gases away from the droplet. Unlike natural 
and forced convection that can be neglected in a spherical case, radial 
convection is always present, although its effects are small at low 
evaporation rates. 

Stefan at 1881 (cited by Fuchs, 1959) was the first to note the 
importance of radial convection. To maintain full pressure in the medium 
together with the partial vapor pressure gradient of the vapor, there must be 
an equal and opposite partial gas pressure gradient of the remaining 
components of the medium. Owing to the presence of the second gradient the 
gas diffuses to the droplet surface, but because of the impermeability of the 
latter the total gas flux towards the surface should equal zero. Hence, the 
hydrodynamic flow of the medium compensates the diffusion of the gas. From 
this discussion it follows that the rate of this flow m* is governed by the 
equation: 
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4.36 


AdC 1 1.1 
o — = m C 
dr 


where C is the concentration and 6 
surrounding gas. Since 


1 


= 6 the diffusion coefficient of the 


C + C 1 = 1 


4.37 


or 


dC 

dr 


1 


dC 

dr 


4.38 


and the total flux of vapor has diffusion and hydrodynamic components, 
Equation 4.27 is replaced by 


dC 


7 r 1 

I = 4Ilp r ( 6 C m ) 

o 1 dr r 


4.39 


dC 

■ 4 <T3J 3^) 


4.40 


Just as Equation 4.27 leads to Maxwell’s equation, so Equation 4.40 gives: 


I q = 4np^ a 6 fcn (1+B); 


4.41 


where the transfer number B is given by 


B = 


C L- C 
1 - C T 


4.42 


The evaporation rate per droplet volume (ra) is given by 


. , 12 P « 

m = I /- Ila = — = £n (l+B) 

o 3 ,2 


4.43 
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4.3.3 Quasi-Stationary Evaporation of Droplets Moving Relative to the 
Media 

The greatest practical interest centers on the evaporation of droplets 
moving relative to the medium under the influence of gravity, inertia, etc. 
This problem can be reduced to the calculation of the rate of evaporation from 
a spherical droplet ventilated by a gas stream. 

Following the majority of the workers in the field of evaporation of 
droplets moving relative to the medium, the convection effects on the 
evaporation rate can be expressed by 

I = I f = I Sh 4.44 

o o 

where Sh is the Sherwood number and f is the wind coefficient. This denotes 
the increase in the rate of evaporation due to the relative movement of the 
medium. 

Using Equations 4.43 and 4.44 one can write: 

1 2<$p 

ra = — ” L *n (1+B) Sh 4.45 

d 

Using the principle of dimensionless analysis one can show that Sh is a 
function of Re and Sc. At the same time the majority of the theoretical and 
experimental work in the field of evaporation of droplets expressed that 
function by 

-nY 

i 

Sh = 2 + B 1 Re i * ' /2 Sc 1 ^ 3 4.46 

where the Sherwood number is given by 
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4.47 


Sh = m /n d 6(C -C) 

Li 

and the Schmidt number is given by 

Sc = v j /6 4.48 

The first accurate measurements of the rate of evaporation of droplets 
suspended in a stream were those of Frossling 1938 (cited by Fuchs, 1959). 

This extremely careful work has served as a model for all subsequent work in 
this field. The experiments were carried out at 20° using droplets of water, 
aniline and nitrobenzene, and spheres of napthalene with a = 0.1 - 0.9 ram 
suspended from glass fibers of radius 25 pm. The droplets were placed 20 cm 
above the exit of a vertical aerodynamic tube of 10 or 20 cm in diameter. The 
Reynolds number was varied over the range of 2.3 - 1280 and the Schmidt number 
range was 0.7 - 2.7. The rate of evaporation was determined by periodically 
photographing the droplet with sevenfold magnification. The temperature 
fluctuation of the air stream did not exceed ± 1%. The determination of the 
rates of evaporation of droplets of organic substances in still air was 
carried out in a closed cylindrical vessel 5 cm in diameter, the walls being 
covered with active charcoal. Frossling examined the effects of different 
factors that might affect the accuracy of the results such as imperfectly 
spherical droplets, turbulence, compressibility of the air, nonideality of the 
vapor, and the counter flow effects and indicated that they are less than the 
experimental error (< ± 1%). His experiments confirm the accuracy of Equation 
4.46 with an experimental value = 0.552. 

The fundamental experiments of Ranz and Marshall (1952) were done at room 
temperature on benzene and water droplets with a = 0.06 - 0.11 cm. The 
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droplets that suspended from the capillary end of a microburette of radius 30 
- 50 Vm were ventilated by dry air from below. The study was restricted to a 
Reynolds number range of 0 to 200. Evaporation rates were determined by 
measuring the rate of feed through a burette necessary to maintain a constant 
droplet diameter. Ranz and Mar hall’s results also confirm the accuracy of 
Equation 4.46 with an experimental value of 3^ = 0.6. 

Ahmadzadeh and Harker (1974) summarized the previous experimental data on 
the evaporation from liquid droplets. All the experimental data in the range 
of Re < 1000 give the value of 3^ =0.55-0.6. 

Kinzer and Gunn (1951) used water droplets of a = 5 - 70 p at 0 - 40°C 
and 10 - 100% relative humidity and employed instantaneous photography to 
measure the evaporation rate. The droplets fell freely in a tube of square 
cross-section. For droplets of a = 0.02 - 0.5 mm the experiments were 
conducted in tubes 200 m long. Insulated metal rings were placed horizontally 

along the axis of the tube at equal intervals. The droplets emerged from a 

\ 

capillary connected to the terminal of a battery and became charged when they 
broke away. Kinzer and Gunn measured the rate of fall of the droplets by 
passing them through a ring, thus creating an electric impulse that is 
amplified and recorded on a moving tape. Since the decrease of terminal 
velocity with time fall gave the rate of mass loss, they determined the rate 
of evaporation as a function of the rate of the droplet fall. They also 
measured the rate of fall and the evaporation rate by a photographic 
technique. For Re = 100 - 1600 their results were fitted by Equation 4.46 
with Sj = 0.46. When Re < 0.9, the constant 3j rises from zero to a maximum 
value of approximately 0.92 at Re = 4 and then gradually falls to the 
abovementioned value of 0.46 at Re = 100. This result contradicts the data of 
other workers, where 3^ has a value in the Re region of 1-100 equal to or less 
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than the value at greater Re. 

Galloway and Sage (1967) have reviewed the available Information 
concerning the effects of the molecular properties of the fluid, conditions of 
flow, and level of turbulence on the evaporation rate from spheres. At normal 
conditions and an intermediate range of Reynolds number. Equation A. 46 
with 0^ = 0.6 represents the standard curve for all other date. 

Yuen and Chen (1978) investigated the evaporation of liquid droplets at 
high temperatures. The data on water and methanol droplets (porous spheres) 
evaporating into the flow of air with the temperature of the latter within 150 
- 960°C was obtained in an air tunnel. The experiments were carried out for 
the Reynolds number range from 200 to 2000. 

Yen and Chen pointed out that for low temperatures their results are 
identical with the standard curve of Equation 4.46 with the same coefficient 
as that of Frossling. 

Prakash and Sirignano (1980) studied the liquid droplet vaporization in a 
hot convective gaseous environment. They developed a new gas-phase viscous, 
thermal, and species concentration boundary layer analysis using an integral 
approach. The gas-phase analysis was coupled with a 1 iquid-phase analysis for 
the internal motion and heat transfer. The coupled equations were solved for 
different hydrocarbon fuels in air at 1000°K and 10 atm. They concluded that 
the heat flux into the liquid phase should be considered in such analysis. 

The temperature distribution inside the droplet is nonuniform for most of the 
droplet lifetime. The Ranz-Marshall correlation seems to agree well when the 
heat flux into the liquid phase is taken into account by modifying the heat of 
vaporization. If the droplet has the same temperature as the surrounding gas 
at the droplet generator exit, the heat flux into the liquid phase can be 
neglected, and in this case the Ranz-Marshall correlation does not need any 
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modifications . 


Now after the discussion of the previous experimental work on the 
evaporation of liquid droplets suspended in a moving stream, it is clear that 
Equation 4.46 as: 

Sh = 2 + 0.55 Re 1/2 Sc 1/3 4.49 

This equation along with Equation 4.45 will be used in the turbulence model to 
calculate the evaporated mass from the droplets to the surrounding gas. 

4.4 Drag Coefficient 

The drag coefficient of spherical solid particles, nonevaporating 
droplets, and evaporating droplets is discussed in this section. 

4.4.1 Drag Coefficient of a Solid Particle 

All the solutions with low inertia terms, Stokes and Oseen, are valid for 
very low Reynolds numbers. The Stokes solution, which ignores completely the 
inertia terms, is valid for Re < 1 ; Oseen considered the Navior-Stokes 
equation with very limited inertia terms, but the drag coefficient is 
unchanged. For a moderate Reynolds number 1 < Re < 200, there are a lot of 
experimental results for the drag coefficient and the plot of these data 
versus the Reynolds number is called the standard drag curve. The recommended 
drag coefficient on a solid sphere in steady motion as the best approximation 
for this curve is given by Clift et al., 1978 

C Dg = (24/Re) [1 + 0.135 (Re)°* 82_0 * 05w J 0.01 < Re <20 4.50 
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4.51 


C DS = < 24/Re) U + 0*1935 (Re) 0,6305 J 20 < Re < 200 

where w = Log^Q Re and the particle Reynolds number is calculated from 

Re = P x |U - Vf d/Uj 4.52 

♦..-22 

u=/u z +u 4.53 

z r 

^=/V 2 +V 2 4.54 

z r 

4.4.2 Drag Coefficient of a Nonevaporating Droplet 

To satisfy the continuity of tangential shear stress across the liquid- 
gas interface, a slight amount of internal motion seems certain to develop. 
This internal circulation of the liquid droplet decreases the boundary layer 
thickness of the exterior flow and may reduce the drag coefficient. 

Pruppacher and Beard (1970) studied the internal circulation and shape of 
a water droplet, which can be considered as a nonevaporating droplet falling 
at terminal velocity in air at 20°C at sea level pressure, and nearly water 
saturated. They concluded that the maximum surface velocity, at the equator, 
of a droplet was found to be about 1/100 of the droplet's terminal velocity. 
Due to this vanishingly small value, one can expect that the flow structure 
around a nonevaporating liquid droplet falling in air will be approximately 
the same as that for a solid particle, hence the drag coefficient will be the 
same. Beard and Pruppacher (1969) measured the drag on water droplets falling 
in water saturated air at terminal velocity in a wind tunnel for Reynolds 
numbers between 0.2 and 200. They concluded that for this Reynolds number 
range the drag on water droplets is in good agreement with that for the drag 
on solid spheres measured or calculated by many other researchers. Beard 
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(1976), depending upon all the available theoretical and experimental data, 
concluded that droplets with diameters < 1 mm (Re < 300) are essentially 
spherical and the drag may be closely approximated by the drag on a rigid 
sphere. Ingebo (1956) investigated the drag coefficient for liquid droplets 
and solid spheres accelerating in air stream using a high-speed camera 
technique. Accelerations of the order of 20,000 m/s^ were considered. The 
sphere diameter range was from 20 to 120 microns. To ensure the spherical 
shape for the liquid droplets (isooctane, water, and trichloroethylene), the 
Reynolds number was in the range 6 < Re < 400. The main purpose of Ingebo' s 
work was the study of the effects of the rate of acceleration, the liquid 
status and the evaporation rate on the drag coefficient. His main conclusion 
was that the unsteady-state drag coefficients are different than the steady- 
state values, but when the acceleration rates were low, the unsteady-state 
drag coefficients were in agreement with steady-state values of previous 
investigations. The interesting result is Ingebo's conclusion that the drag 
coefficient for slowly evaporating droplets, nonevaporating droplets, and 
solid spheres are the same. 

Rivkind et al . (1976) solved the Navier-Stokes equations for the flow of 
fluid inside and outside a spherical droplet using the method of finite 
differences. They considered the variables Re, ^ anc * P 2^ P 1 aS contr °m n 6 
parameters of the problem. The drag coefficients of the droplets were 
investigated for 0 < ^^1 < 00 and 0- 5 * Re * 100. They concluded that the 
density ratio has virtually no effect on flow characteristics. According to 
their numerical results, they proposed that the drag coefficient of the 
droplet can be defined in terms of the drag coefficients of the solid sphere 
(C DS ) and of the gas bubble (C Db ) at the same Re by the formula 
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Rivkind and Ryskin (1976) extended their previous work to consider the 
circulating flow inside and outside the droplet up to Re = 200. They again 
recommended Equation 4.55 for the calculation of the drag coefficient of a 
liquid droplet moving in gaseous flow. By calculating the drag coefficient of 
water and methanol droplets moving in atmospheric air with a Reynolds number 
up to 200, it was found that the difference between the values produced by 
Equation 4.55 and those of a solid particle (Clift et al . , 1978) is only 3%. 

Hamielec and Johnson (1962) used the error distribution method to predict 
the velocity profiles and terminal velocities for solid and fluid spheres 
moving in viscous media under the influence of gravity. The error 
distribution or Galerkin method is based on choosing a polynomial for the 
stream function that satisfies all the boundary conditions together with an 
integral form of the Navier-Stokes equation. By this method, Hamielec and 
Johnson predicted reasonably accurate velocity profiles and terminal 
velocities for circulating droplets and bubbles. Hamielec et al. (1963) 
modified the work of Hamielec and Johnson (1962) to account for the finite 
interfacial interface with trial stream functions. They predicted velocity 
profiles for viscous, laminar, and imcompressible flow around droplets, 
bubbles and solid spheres. The drag coefficients, flow separation angles, and 
forced convection transfer rates were calculated and compared with 
experimental data for solid spheres, circulating droplets, and bubbles of some 
other workers. They obtained good agreement up to Re = 100. Hamielec et al. 
(1963) tried to correlate the available experimental data for the drag 
coefficient of liquid droplets falling in water using a viscosity-ratio 
correction factor. This relation is given by: 


4.56 


„ 3.05 , 78 V + 214 2Ur+ 1080 , 

D ~ Re?* 74 " 1(60 + 29 u r ) (4 + 31^) j 
10 < Re < 100 

where 

In the present study the drag coefficient of a water droplet moving in 
atmospheric air was compared with a Reynolds number in the range 
(10 < Re < 100) using Equation 4.55 and 4.56 with the experimental data of 
Beard and Pruppacher (1969). It was found that the values produced by 
Equation 4.55 are 3% less than the corresponding experimental values while 
those produced by Equation 4.56 are 17% less. Therefore, Equation 4.55 is 
recommended for the calculation of the drag coefficient of a moving 
nonevaporating droplet in a media with comparable viscosity to that of the 
droplet. 

Nakano and Tien (1967) also used Galerkin's method to investigate the 
effect of increasing the internal Reynolds number (0 < Re < 50) or, more 
accurately, the flow behavior within the fluid sphere by including inertia 
terms for both phases. Changes in the internal Reynolds number had little 
effect on the external streamlines or on overall drag. Thus, they got almost 
the same results as Hamielec and his coworkers. 

Now, based on the previous theoretical and experimental work, it is clear 
that the effect of internal circulation and low evaporation rate on the drag 
coefficient of a liquid droplet moving in a gaseous stream is negligible. 

This drag coefficient can be considered from the standard drag curve of a flow 
over a solid sphere. 

4.4.3 Drag Coefficient of an Evaporating Droplet 

Hamielec et al. (1963) studied numerically the effect of mass transfer 
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from a spherical particle at Reynolds numbers 1, 40, and 100 on the drag 
coefficient. In all cases they showed that radial mass efflux decreases 
friction drag and increases the pressure drag slightly. Due to the bulk flow 
of vapor from the droplet surface, the boundary layer thickness decreases. 
Thus, the velocity gradient or the surface shear stress decreases, so a 
reduction in the friction drag is predicted. The slight increase in the 
pressure drag may be attributed to the blowing effects on the angle of 
boundary layer separation. On the other hand, Kassoy et al. (1966) have shown 
that at low Reynolds number, the drag of a sphere at constant free stream 
temperature decreases with decrease in temperature of the sphere. This is 
attributed to the changes in the surrounding fluid properties that might be 
more pronounced in the case of evaporation. Yuen and Chen (1976) have noted 
that the changes in the composition of the gases near the droplet surface are 
important and would tend to reduce the drag of an evaporating droplet since 
the viscosity of most vapors is lower than the viscosity of air at the same 
temperature. The effects of the temperature and concentration gradients due 
to evaporation on the dependence of drag coefficient on Reynolds number are 
accounted for by using the free stream density and the 1/3 rule for the 
dynamic viscosity (Yuen and Chen, 1976). 

Eisenklam et al . (1967) investigated experimentally the evaporation rates 
and drag coefficients for evaporating and burning droplets of various fuels 
freely falling in atmospheric air at temperatures of up to 1000°C, or burning 
in cold oxygen atmospheres. They correlated the experimental data for the 
drag coefficient Cp using the formula 


C D = 


C DS /(1 


+ B) 


4.57 
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The experimental data covers the Reynolds number range 0.005 - 15. The 
correlation 4.57 was suggested by results from boundary layer theory and an 
analogy with droplet heat transfer, for which case a similar correlation was 
already in existence. 

For intense mass transfer, such as when the droplet is burning, the 
evaporation is expected to reduce film drag due to a thickening of the 
boundary layer, and if the droplet is burning, the form drag is expected to be 
reduced by the "filling in" of the wake by products of combustion. Further 
effects are the alteration of the position of boundary layer separation and 
the steep variations in properties due to the large temperature and 
concentration gradients associated with intense mass transfer. 

Yuen and Chen (1976) used Eisenklam's experimental data, along with their 
own data that extended the Reynolds number range to about 500, to develop an 
alternative correlation. They defined a reference Reynolds number 

= Re « y «/ u r 4 * 58 

where is the viscosity of a reference mixture at the temperature 

T = T t + (T - T t )/3 4.59 

r L 00 L 

N» 

and containing a vapor mole fraction (Xj.) given by S 

X r -X L +(X a) -X L )/3 4.60 

The subscripts L and 00 denote, respectively, the conditions at the droplet's 
surface and far away from it. Their correlation for the drag coefficient is 
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4.61 


then given by the rule that states that if 
C DS * F <*0 

is the standard drag coefficient curve for solid spheres as a function of free 
stream Reynolds number, then 


C D = F(Re r ) 4.62 

will be the drag coefficient of an evaporating spherical droplet. 

Dukowicz (1984) calculated numerically the drag coefficient of 
evaporating droplets in the Stokes flow regime. He tested Yuen and Chen 
(Equation 4.62) and Eisenklam (Equation 4.57) correlations in the low Reynolds 
number range. Dukowicz (1984) and Sirignano (1983) recommended the use of 
Yuen and Chen's correlation for the calculation of the drag coefficient of an 
evaporating spherical droplet. 



4.5 Effect of Free Stream Turbulence on Drag and Evaporation Rate 
The flow conditions about the particle, especially the free stream 
turbulence intensities, may cause large variations in the drag coefficients 
from those values given by the standard drag curve. Zarin and Nicholls (1971) 
reported that at Re < 200, they have observed little or no change in Cpg 
compared to the C D g's measured at lower turbulence intensities ( « 1%). 

In general, the motion of a droplet in a turbulent flow depends upon the 
characteristics of the droplet and of the turbulent flow. Droplets with small 


size or small realization time ( t ) compared with the turbulence time 
scale ( t ) respond to the fluctuating motion of the carrier fluid. If 

L 

t > , very little fluctuation in the droplet velocity can be seen (Clift 
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et al . , 1978). In this case the effect of turbulence is then to modify the 
flow field around the droplet, so that the drag may be affected. 

In the present study the value of the ratio x / x T is expected to be less 

p L 

than unity, and thus the effect of the free stream turbulence on the drag 
coefficient should be very small. Experimental evidence supports this 
assumption. For example, Clift et al. (1978) reported that the effect of the 
free stream turbulence in the range of droplet Reynolds number 10 < Re < 50 
is less than 5% (Figure 10.11 in that reference). 

Regarding the evaporation rates, for Re < 50 and for turbulence intensity 
less than 20%, the experimental data (Clift et al . , 1978) showed a very small 
increase in the evaporation rate. Therefore, in the present study the drag 
coefficient and the evaporation rate will be mainly functions of the Reynolds 
number only and not in terms of the turbulence intensity around the particle. 
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5.0 EDDY DIFFUSIVITY OF A SINGLE PARTICLE 

5.1 Introduction 

Eulerian mathematical models of particle-laden turbulent flows require 
knowledge of the statistics of particle motion in order to calculate the 
particle response to fluid turbulence. Specifically, a reliable expression is 
needed for calculating the Schmidt number, defined as the ratio of particle 
diffusivity to fluid point diffusivity for heavy particles. So far, such an 
expression is unavailable in the literature. 

Tchen (1947) studied the diffusion of a rigid particle carried by a 
turbulent flow. Hinze (1975) indicated that Tchen's assumption that the fluid 
element should continue to contain the same discrete particle at any time is 
hardly to be satisfied. Such an assumption cannot be valid if the ratio 
between the material densities of the particle and fluid is large. In this 
case the particle will be associated with more than one eddy along its path 
which is termed the overshooting phenomena. Soo (1956), Friedlander (1957), 
Chao (1964), and Gouesbet et al. (1984) solved Tchen's equation under the 
assumption of no overshooting. For asymptotic times of dispersion, that 
assumption implies a particle Schmidt number of unity which is physically 
incorrect . 

Peskin (1971) studied the particle diffusivity under the condition of 
overshooting but restricted his analysis to small distances between the 
discrete particle and the "originally surrounding fluid." As a result, 

Peskin' s formula predicts values of particle Schmidt number very close to 
unity. Some other workers (Reeks, 1977, and Nir and Pismen, 1978) 
investigated the particle diffusivity under the condition of overshooting but 
restricted their study to the Stokes flow regime, which is hardly to be 
satisfied for suspended heavy particles in a gaseous media. 
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Csanady (1963) studied the differences between the diffusion of the fluid 
points and heavy solid particles in the atmosphere. He attributed the 
appreciable reduction of the dispersion rates of the heavy solid particles to 
the rapid travel across the turbulence eddies. Meek and Jones (1973) 
statistically analyzed the heavy particle behavior for a constant relative 
velocity between the particle and its surrounding gas in homogeneous 
turbulence. They obtained expressions for the particle dispersion coefficient 
and its Lagrangian autocorrelation. 

In this section a reliable expression for the calculation of the lateral 
diffusivity of heavy particles suspended in a homogeneous turbulent field is 
provided. The physical parameters that control the particle behavior in a 
turbulent flow is discussed. Two widely used theories (Csanady, 1963, and 
Meek and Jones, 1973) for the calculation of the statistical properties of 
heavy particles in terms of those of the surrounding fluid are reviewed. The 
empirical coefficient in Csanady* s theory is determined via a comparison with 
the experimental data. Finally, Meek and Jones' theory is examined and an 
empirical coefficient is introduced and evaluated via a comparison with the 
experimental data. 

5.2 Physics of Particle Dispersion 

The behavior of a single spherical particle suspended in a homogeneous 
isotropic turbulent field depends on the properties of both the particle and 
the turbulent flow. The first parameter that controls the particle dispersion 
is the ratio of the particle size to the Kolmogorov, length scale, n. If this 
ratio is small the particle dispersion will be influenced by the entire 
spectrum of eddy sizes and will follow the turbulence fluctuations of the 
carrier fluid. If, on the other hand, that ratio is greater than unity, the 
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particle will follow the slow large-scale turbulent motions of the fluid 
(Alonso, 1981). In this case the main effect of the turbulence on the 
particle is to modify the flow field around it, so that its drag may be 
affected. 

A more general parameter is the ratio between the particle relaxation 

time, x , inversely proportional to the particle's inertia, and the fluid 

Lagrangian integral time scale, t , is the controlling parameter of the 

L 

particle response to the turbulence fluctuations. If the ratio t /t < 1, the 

pi* 

particle will be able to respond to the entire spectrum of fluid motion and 

the ratio between the root mean square fluctuating velocity of the particle to 

that of a fluid point is almost unity. On the other hand if x /t > 1 the 

p L 

particle will respond very slowly to the fluctuating fluid motion. 

If there is an appreciable relative mean velocity between the particle 

and the surrounding fluid, the particle will move about from eddy to eddy, 

. L > \ 

whereas a fluid point would remain in the same eddy throughout the lifetime of 
that eddy. This is what Yudine (1959) called the "crossing trajectories 
effect." Yudine (1959) considered the physical consequences for finite free- 
fall velocities on the heavy-particle diffusion. He formulated upper and 
lower limits for the changes in the dispersion coefficient due to the heavy 
particles free fall velocity, f. 

Taylor (1921) postulated a theory describing the statistical dispersion 

of fluid points in a stationary homogeneous turbulent flow. His result 

~2 

relates the mean lateral square fluid point displacement, Y , to the mean 

~~2 

square fluctuating velocity, u , and to the Lagrangian velocity correlation 
coefficient, R L (t), according to: 

Y 2 (t) = 2 u 2 0 J t Q J T R l (t) dtdT 5.1 
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The correlation coefficient is defined by 

R (T) = u(t) u(t + T) 

L ~1 
u 

Taylor (1921) defined a turbulent diffusivity by 
E t (t) -rjt y2 (t) 

= u 2 0 j T Rjt) dT 5.4 

■ “ 2 \ ■ 5 - 5 

Snyder and Luraley (1971) showed that Taylor's theory is equally 
applicable for the dispersion of alien particles, provided that the velocities 
are interpreted as particle velocities, and that the Lagrangian time 
scale, T^, is interpreted accordingly. Snyder and Lumley (1971) and Wells and 
Stock (1983) measured the dispersion of heavy particles in a grid-generated 
turbulent flow. They reported the decrease of the particle integral time- 
scale and the rapid decrease of the particle velocity correlation with 
increasing t . 

Calabrese and Middleman (1979) photographically measured the degree of 
radial dispersion of medium-size particles emanating from a point source in 
the turbulent core of a fully developed vertical pipe flow of water. They 
defined the medium-size particles by h < d < , where d and are the 

particle diameter and Lagrangian length scale respectively. They were able to 
calculate the radial mean-square particle displacement directly to find that 
both heavy and buoyant particles experienced a decrease in mean-square 


5.2 
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displacement due to the crossing-trajectory effects. 


5.3 Csanady's Theory 

Csanady (1963) considered the differences between the diffusion of the 
fluid points and heavy particles in the atmosphere. He proposed a functional 
form for u(o) u(t), where u(t) is the fluctuating velocity of the air, 
consistent with similar shapes for Eulerian and Lagrangian fluid point 
correlations. Thus, Csanady was able to construct two relations for the 
particle dispersion coefficient parallel to and normal to the direction of the 
particle terminal velocity. Lumley (1978a) showed that Csanady's model gives 




1 


(1 


♦ s 2 fV )" 2 


6 h - 4/3 


5.6 



(1 + B 2 f 2 /u 2 ) 1/2 

V 


% - 2/3 


5.7 


where is the asymptotic dif fusivity ' for the fluid. The subscripts h and v 
correspond respectively to horizontal and vertical dispersion of the particle 
associated with long times relative to the integral time scale of the 
turbulence. Lumley stated that the value of ^ is much less well determined 
than that of 8^ and the value of 8^ should be determined via a comparison with 
well defined experimental data. 

In the present work the value of 8^ is determined by comparing the 
prediction using Equation 5.6 with the experimental data. 
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5.4 Meek and Jones' Theory 

Meek and Jones (1973) considered the motion of a heavy particle in a 
homogeneous turbulent air flow. The particle motion was viewed from a 
reference frame in which the average fluid velocity is zero. They started 
with the definition of the particle velocity autocorrelation, which can be 
expressed in terms of a normalized particle energy spectrum, E_ ( w) : 

R ,,( T ) = lZ E ,.(<*)) cos (ojt) doi 5.8 

where w is the circular frequency. 

Soo (1956) obtained an expression relating E n jj (id ) to its fluid 

P > x O 

counterpart, E ii ( a, )> from the solution of the simplified particles equation of 
motion and under the assumption of zero free-fall. This expression is given 
by 


E p,U ( “o ) ’ W ‘ll^' 


5.9 


where is the particle response function defined as: 


u 


Q ii ( V = 


ii 


ii 


1 + (-—) 
°1 


-] , 


5.10 


18 

where <*. = ; 

V 


5.11 


d and P2 are the particle diameter and density respectively and is the 
dynamic viscosity of the surrounding fluid. The Lagrangian frequency 
function, (u>), can be approximated by various semiempirical forms. Hinze 
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(1975) pointed out that the use of an exponential form for the fluid's 
Lagrangian velocity autocorrelation provides reasonable agreement with the 
experimental data. The associated energy spectrum is 


E ii<"> ■ 


L,ii 


' O * J'ltP 


5.12 


Meek and Jones (1973) indicated that the use of a response function 
derived for zero free-fall velocity, f^, requires some adjustment of 5.8 to 
account for non-zero free-fall and the subsequent movement of the particle 
from one eddy to another. They suggested that the nonezero free-fall velocity 
spectrum E p ^ ( u>) , be expressed in terms of the zero free-fall velocity 
spectrum, according to: 


E n <;(■«>■ = E 5 

P.ii P.ii o 


5.13 


where Z = [ 1 + (f^/v^)^]^ 


5.14 


Using the above assumptions, the particle velocity autocorrelation can be 
written as 


R „(■*) - 

p ,ii 


1 - 5 


[e T/T li - 


ii 


5 it i T/C ii T ii], 


5.15 


where 5^ 


1 L,ii 


5.16 


T ii = X L,ii /Z * 


5.17 
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The ratio of the fluctuating velocity variances is given by Soo (1956) as: 


(— ) = 


1 


2 

U 


5.18 


while the mean square particle displacement, Y^(t), in the radial direction is 


~2 v2 T n 

Y^( t) = 11 

P 


y _ | — [(1 - e t/T ll) - ( 1 - i t/(5 ll T ll ) )] . 5.19 


’ll 


The particle eddy diffusivity is given by 


V 2 T 

e ( t ) = [(i _ e t/T ll) - 5 2 (1 - e t/(5 ll T ll ) )] . 

% 1’1 


5.20 


The corresponding fluid point eddy diffusivity is given by 


e f (t) = u 2 t l>11 (1 - e t/T ll) 


5.21 


Thus, the particle Schmidt number (a = — ) is given by: 

p e 


e h (t) 

°p (t) " W 


d - 5 u )"z 


r , t/T lK 2 , ^1 1 T ll ) 

[(1 - e ) - K ■ i (l - e 


] 5.22 


For long times relative to the integral time scale of turbulence 


a = J 1=1, 

P e f Z ’ 


5.23 


where Z is given by (5.13). 
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Peskin (1974) argued that a heavy particle initially coincident with the 
fluid point naturally lags behind that fluid point as a result of inertia. 

The fluid points encountered by the heavy particle are not statistically 
equivalent and would only be so in a turbulent field with infinitely large 
space and time Eulerian correlation scales. Finite correlation scales imply 
that the fluid point encountered by the heavy particle as it lags behind the 
original coincident point are not statistically identical. Accordingly, the 
heavy particle dispersion or diffusivity is determined by the Eulerian space- 
time correlation. Peskin stated that Meek and Jones failed to account for the 
most important effects of the space-time correlation of the turbulent flow on 
the particle motion. The linear relation (Equation 5.9) between the particle 
energy spectrum and the fluid Lagrangian spectrum implies that the effects of 
the space-time correlation are completely neglected. 

Meek and Jones (1974) in their reply to Peskin (1974) argued that the 
good agreement between their predictions and the experimental data of Snyder 
and Lumley supports their solution, especially in homogeneous turbulence. To 
achieve this agreement, they used Equation 5.19 and the fluid and particle 
data of Snyder and Lumley (Equation 5.28 and Table 1 in Meek and Jones' 
article, 1973). Using their equation and table, the present study can not 
reproduce their Figure 2, especially for the solid glass particles. 

5.5 Modifications of Meek and Jones' Theory: 

As pointed out by Peskin (1974), the Eulerian space-time correla- 
tion Rjj,(y(x) - x( t) ) , where y( x) and x( x) are the Lagrangian fluid and the 
particle positions, should be considered in the analysis of heavy particle 
dispersion instead of the Lagrangian autocorrelation. To do so the solution 
should be restricted to a very short distance (y( x) - x(x)) between the 
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particle and the fluid, otherwise the equations will be formidable. Peskin 
imposed this restriction but his solution predicts values of a particle's 
Schmidt number very close to unity even for a heavy particle. Yudine (1959) 
in his discussion of the physical parameters controlling the heavy particle 
diffusion pointed out that the dispersion process is controlled mainly by the 
terminal velocity, f. He stated that the dispersion process depends upon the 
terminal velocity in three ways: (1) the terminal velocity determines the 

vertical displacement of the center of dispersion of the particle; (2) because 
the terminal velocity is a certain measure of inertia, the particle does not 
follow completely the high frequency fluctuations of turbulent fluid velocity; 
and (3) if it has appreciable terminal velocity, a particle will fall from one 
eddy to another, whereas a fluid' point would remain in the same eddy 
throughout the lifetime of the eddy. Yudine concluded that, for large f, the 
dispersion coefficient takes on an asymptotic form inversely proporatlonal to 
f. 

Meek and Jones (1973) pointed out that the inertial effects can be 
significant especially when f^ < v-^. The inertial effects increase the 
particle Lagrangian time scale compared to that of the fluid if there is no 
crossing trajectories. 

Csanady (1963) accounted for the crossing trajectory effect on a heavy 
particle dispersion by including in his analysis the quantity f^/u^ and the 
ratio between the Lagrangian to Eulerian integral time scales, 6^. 

Due to the close similarity between the two theories of Csanady and Meek 
and Jones and based upon the previous discussion it is clear that the 
parameters that should be considered in a heavy particle diffusion analysis 
are the terminal velocity and a coefficient simulating the ratio between the 
Lagrangian to Eulerian integral time scales. Therefore the following 
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modification to the Z factor (Equation 5.14) is proposed: 


1/2 


Z - a + V 2 i/v j> 


5.24 


where <x^ is an empirical coefficient to be determined via comparison with the 
experimental data. 

To determine 8^ in Equation 5.6 and in Equation 5.24 the predictions 
will be compared using both Csanady's theory and the modified theory of Meek 
and Jones with the experimental data of Synder and Lumley (1971) and Wells and 
Stock (1983). Since the two theories have been developed for stationary, 
homogeneous turbulence, the data should satisfy these two conditions. 
Experimental evidence suggests that grid-generated turbulent flows approximate 
the stationary requirements and corrections can be made for the inhomogeneity 
(Pismen and Nir, 1979). 


5.5.1 Snyder and Lumley' s Experiment 

Snyder and Lumley (1971) performed an experiment in which a single 

spherical solid particle was injected into a turbulent flow generated by a 

grid. They considered particles of various sizes and densities ranging from 

light particles that closely follow the fluid fluctuations to heavy particles 

that experienced both inertia and crossing-trajectory effects. The particles 

were injected above the gird at a distance of 20 grid spacing. They measured 

the particles' Lagrangian autocorrelations at x/M = 73 as well as the mean 

~2 

square displacement, , of the particles as they were individually convected 
through the wind tunnel . Turbulence measurements were made and the turbulence 
energy decay was given by 
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5.25 


U 2 x 

— = 42.4 (£ - 16) 

T M 

u 

where U is the mean axial velocity (6.55 m/s) and x/M is the ratio of the 
axial distance to the mesh size (2.54 cm). They also reported the rate of the 
energy dissipation, £, along the centerline of the wind tunnel (see Table 
5.1). 


Table 5.1. Energy Dissipation Rate ( e) of Snyder and Lumley (1971) 


x/M 

41 

64 

73 

107 

138 

171 

„ 2, 3 

£ cm /s 

5430 

1610 

1160 

480 

266 

165 


Table 5.2 lists the relevant characteristics for the particles studied by 
Snyder and Lumley. 

Table 5.2. Particle characteristics for the 
data of Snyder and Lumley (1971) 


hollow glass corn pollen glass 

3 

density, p g/cm 0.26 1.0 2.5 

2 

diameter, d vm 46.5 87.0 87.0 

drift velocity, f cra/s 1.67 19.8 44.2 


Snyder and Lumley estimated the Lagrangian time scale, t , at one station 

JJ 

only (x/M = 73) by considering the light particle results as representative of 
Lagrangian correlations. In comparison with Snyder and Lumley's data, the 
variation of t along the wind tunnel axis should be estimated from the 

Lt 

measured turbulence quantities. x can be obtained from the relation (Hinze, 

lj 

1975) 
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H-c“ 2/ ‘ 5 ' 26 

The ratio between the Lagrangian and Eulerian. time scales, as given in 
Equation 5.26, depends on the local flow conditions and cannot be considered 
as constant. Hinze (1975) pointed out that the experimental data show that 
the value of the constant C varies between 3 and 10 depending on the Reynolds 
number, thus (5.26) reads 

t l =C’u 2 / £ , 5.27 

where C' varies from 0.2 to 0.66. 

On the other hand, Calabrese and Middleman (1979), using Taylor's theory, 
obtained the following expression: 

x L - if - 0.625 u 2 / e 5.28 

Berleraont et al. (1982) considered the closure relations of the K-g model and 
obtained 

x T = 0.2 u 2 / e 5.29 

L 


In the present work the coefficient C' in Equation 5.27 will be 

determined from the dispersion data of the light particle (hollow glass) or 

the particles with zero terminal velocity that are considered as 

representative of Lagrangian correlations. Using Equations 5.18, 5.19 and 

~2 

5.26 for the hollow glass particles, is obtained and plotted in Figure 5.1 
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with an optimized value of C' = 0.35. Figure 5-1 also compares the prediction 
using Meek and Jones' theory and its modified form with the experimental data 
for the mean-square displacement of the corn pollen and the glass particles. 

It is clear from this figure that Meek and Jones' theory does not provide good 
agreement with the experimental data. With the present modification of Meek 
and Jones' theory the present study predicts the experimental data for the 
particles with different diameters and densities in Figure 5-1. This 
agreement is obtained with = 0.3. 

Figure 5-2 shows a comparison between the theoretical and the 
experimental data for the Lagrangian autocorrelations for the different 
particles. The correlations decrease faster for the heavier particles (high 
drift velocity compared with the turbulence intensity) due to the crossing 
trajectory effects. It is also clear that the value of =0.3 produces 
good agreement with the data for the autocorrelations. 

5.5.2 Wells and Stock's Experiment 

The effects of "crossing trajectories" and inertia on the dispersion of 
particles suspended in a field of grid-generated turbulence were investigated 
experimentally by Wells and Stock (1983). The flow conditions and grid size 
and shape were very similar to those used by Snyder and Lumley (1971) except 
that the main direction of the flow was horizontal instead of vertical. The 
particles were glass spheres, with a diameter of 5 or 57 pm and a density of 
2.45 gm/cm^. The particles were charged by a corona discharge then injected 
on the centerline of the flow. The test section was subjected to an electric 
field, which provides a coulomb force to the particles to balance the 
gravitational force. In this way, the drift velocity could be changed. The 
particle concentration and velocity were measured with a laser-Doppler 


97 







FIGURE 5-2 LAGRANGIAN AUTOCORRELATIONS 
FOR THE DIFFERENT PARTICLES 
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anemometry system. The data were reduced to yield the particle mean-square 
displacement . 


The measured turbulence intensity, ^ , and dissipation rate of kinetic 
energy, £, were used to calculate t . Equation 5.27 was used with C' = 0.5 

Lt 

to reproduce the experimental data for the particles with f^ = 0 as shown in 

~2 

Figures 5-3 and 5-4. Figure 5-3 shows versus x/M for 5 pm particles with 
two values for the terminal velocity: zero and 25.8 cm/s. The two figures 
show good agreement between the predictions and the experimental data 
using = 0.3 that was optimized for Snyder and Lumley's data. 

Figure 5-3 displays the distribution of the Schmidt number, o^, against 
fi/uf This figure compares the prediction using 5.24 with = 0.3, 
Csanady's model (Equation 5.6), Snyder and Lumley's data, and Wells and 
Stock's data. The diffusivity ratio for the experiments of Snyder and Lumley 
and Wells and Stock were obtained using Equation 5.3 at different times or 
x/M. The agreement between this work and the experimental data is very 
good. The solid line in figure 5-5 was obtained using a value of in 
Equation 5.6, equal to 0.55. 

Comparison of the modified form of the Meek and Jones ' theory (Equation 


5.24) and that of Csanady (Equation 5.6) shows both similarities and 
differences (Mostafa and Elghobashi, 1985b). In Csanady's work, the only two 
parameters controlling the dispersion of a heavy particle in a turbulent flow 
are (f/v) and In that way the ratio f/u is considered as a measure of the 

crossing-trajectories effects, together with the associated continuity 
effect. In the present work, the ratios (f/u) and (v/u), and are the 
controlling parameters. The ratio v/u is a direct measure of the inertia 
effect on the dispersion process as discussed by Meek and Jones and other 
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FIGURE 5-3 PARTICLE DISPERSION, 5fA PARTICLES 
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workers (see Hinze, 1975) and it has a significant effect, especially when 
f/v <1. 


5.5.3 The Final Expression for Particle's Schmidt Number 
The expression for calculating the Schmidt number of a heavy particle 
suspended in a turbulent flow (Equations 5.23 and 5.24) has been developed 
assuming that the drift velocity is constant or is large enough to make the 
Eulerian space-time correlation approximated by the Lagrangian autocorrelation 
(Peskin, 1974, and Reeks, 1977). Therefore the developed expression can be 
used for suspended particles in a turbulent flow if the relative mean velocity 
is assumed to be constant during a period of time t > t^. At dispersion time 
greater than the Lagrangian time scale of turbulence, that expression is given 
by 


a = [1 + 0.3 (U - V ) /(vV] 
P 


k 2 

(^_) = 1 


— k * 

2 1 + T K / T 

u p L 


T = 


P 2 (dV 


18 Pj ’ 


t l = 0.35 u / e 


5.30 

5.31 

5.32 

5.33 


The coefficient 0.35 in Equation 5.33 is the optimized value using the 
experimental data of Snyder and Lumley (1971). 
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6.0 NUMERICAL SOLUTION OF THE EQUATIONS 

We have now available a turbulence model and a number of experimental 
target values for free jet flows. The present section shows how the model 
will reproduce these target values. 

To apply the model, a solution method and boundary conditions are 
needed. The first part of this section disucsses the of finite-difference 
technique used in solving the differential equations. The second part 
disucsses the solution procedure and the boundary conditions. 

6.1 The Equations to be Solved 

At this point, the reader may welcome a brief reminder of what 
constitutes the prediction method applied in this section. 

This method employs the mean equations of continuity and momentum for 

each phase, the global continuity equation, the concentration equation, and 

the K and e equations. Thus, the equations to be solved are (6 + 3 k) in 

k. k 

number; they are for the dependent variables U z , U r , , V r , 

P, ^ K and e * 

6.2 Solution Method 

The governing equations are solved simultaneously with the finite- 
difference method that Spalding (1979) has developed for laminar two- 
dimensional parabolic dispersed-f low problems with interphase slip (GENMIX- 
2P). Since the governing equations are parabolic in nature, the method 
integrates by marching forward, i.e., downstream, starting at an initial 
cross-section where the profiles for all dependent variables must be 
specified. The GENMIX-2P computer code is generated from the GENMIX computer 
code (Spalding, 1978), after excluding the effects of mass transfer, chemical 
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reaction, and turbulence then adding two new subroutines. These two 
subroutines are C0MP2P and ADJ2P, which set up and solve the finite difference 
equtions for the dispersed phase variables. All the necessary information 
about GENMIX is documented (Spalding, 1978) and need not be repeated here. 
Therefore, we will restrict the description here to the treatment of the 
dispersed phase and how it is coupled with the carrier phase in the solution 
procedure. 


6.2.1 The Computational Mesh 

The computer code GENMIX-2P employs the stream function ^ of the carrier 
phase as a cross-stream variable that is defined as 

* - J 0 r pi*iV dr • 6,1 

The governing equations are transformed into a coordinate system based on the 
axial distance, z, and a normalized stream function, .defined as 



where w and y are the values of the stream function at the inner and outer 
I E 

boundaries of the flow, s measures the distance from an arbitrarily assigned 
starting point that is often taken as the starting point of the marching 
integration along the inner boundary of the grid. Lines of constant z are 
normal to the I boundary, is assigned the value 0.0 along the I boundary 

of the grid. It’s value then increases monotonically with distance from that 
boundary, rising to the value 1.0 at the outer boundary of the grid. 

The flow field of interest is subdivided into cells and each cell is 
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treated as a control volume. For the carrier phase, these cells are formed by 
the intersection of the constant z and oj^ lines as shown in Figure 6-1. 


6.2.2 Finite Difference Equations (FDE) of the Dispersed Phase 

The finite-difference equations for the carrier and dispersed phases are 

formed from the differential forms by integration over the control volumes of 

the carrier phase as shown in Figure 6-2. In this figure UP denotes upstream 

station and DN stands for downstream. TE^ is the tangent of the 

constant line just above the grid point i, but it is not distinguished from 

the angle or its sine. oG^ is the inclination of the streamlines of the 

dispersed phase at i ~ i + 1 interface. For simplicity, in this section only, 

one class of particles will be considered. Therefore, the superscript k will 

be replaced by the subscript 2 to represent the dispersed phase. al^ equals 

TE^ plus an increment allowing for cross-flow, which takes account of 

local p. $ U . Therefore, it represents the tangent of the streamline angle, 

1 J. z 

but it is not distinguished from the angle or its sine. is the outlet flow 

f~ h 

area of the i L control volume and Az = z^ - z^. 

For the variables I . and V I . the control volumes over which 

2 ' i z 1 i 

integration of the equations is carried out are those of U z l^ of the carrier 
phase. For the control volume is the one bounded by i + l,i at the 

upstream and downstream stations. Let Q represents any variable of the 
dispersed phase such as V 2 , V r or a2 and ^ • The result of integrating the 
equations can be most conveniently expressed as follows, 


Qi = 


A i Q i+l + B i Q i-l + C i 

D. 


6.3 


The and coefficients represent the effects of transport of the dispersed 
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direction of 
marching integration 





Figure 6-1 ILLUSTRATION OF THE z-W COORDINATE SYSTEM 
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I 



Figure 6-2 ILLUSTRATION OF THE CONTROL VOLUMES USED 
FOR THE DERIVATIONS OF THE FINITE DIFFERENCE EQUATIONS 
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phase across constant o> n lines. The coefficient represents the effects of 
upstream convection from the different sources. These sources result from the 
pressure gradient, the gravity force, the interphase friction in the case of 
V z and V f , and the different turbulent correlations that result from the 
presence of the particles. The coefficient represents the effect of 
outflow from the control volume. In the following the forms of the 
coefficients A^, B^, and for the variables a2 (=V r /V z ), V z and are 
presented: 

For a2: 


=> max [0, (DFa - 0.5 CONcO^^^* “ CONa i+3/2^ V zl i+3/2 * 6.4 

B t = max [0, (DFa + 0.5 CON cx> i+1 / 2 ^ V Z ^ i+l / 2 » 6,5 

C i = p 2 (4 2 SV z V z ) li+l/2 a2 il + (P i~ P i+l )r i+l/2 Az 

*2 I i+ 1 / 2~ g r^ p 2 — P 1 ^ ®2 1 i+1/2 V0L + FU z®2 I i+l/2 otl i V01, » 6,6 

D i = lp 2 $ 2 SV Z li+l/2 + A i /V z I i+3/2 + B i /V z ^ i+1/2 + 

Vli + 1/2 V0L] V z I i+1/2 ’ 6 * 7 

where VOL is the total volume of the cell and, 

CON V 3 /2 ■ <r Vz )| l+l (ce - IE) 1+1 ■ 6 - 8 

DFa i+3/2 ’ - I'V'cO - 1/0 p )4 2,r 1 i + l/2 /V 2 l 1+ l V0L 
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(1 - TE/a2) i+ i 


- Po 62 H+l 1 p 


V'i*‘ 


+ Pjto r i + i*2 l i+i lv + P 


2 + » V z r 1 El 1 + 1 / 1 + 3/2 


- p,Az l i+l 1 P 


r ,., 1 V°C ( * 2 ,r ' 2 - 


- *, te)] 1+1 


+ P 2 te c 4 > 


r ~ 2 — V V TEl 

trK/e (v t /o c » 2jI ) jt 12 P z - r 


1 


• *_ r 0 O / v „ I i + 1 

+ m ® r i+l p c 2 z 


6.9 


For V : 


= max [O.CDFV 


-0.5 C0NV) i+1 / 2 * * c0NV i+l/2 


1 , 


. max 10, (OPV ♦ 0.5 COW^-l/lI ’ 


6.10 


6.11 


C t - ■Pjl S *i V iVi " <Pd ‘ P,,) l * 2Sli 

82 (V»P VOL » 2 U + l F Vx'i 

„ t . ^isv/i + ^ * lFVl V ° L ’ 


6.12 


6.13 


where 


CONVj - P 2 ^ tv¥x' 1.1/2 


(ofi - TE) i 


6.14 
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DFV. - - 4 « (« 2 ) 1+ , /2 lv p O v l I+ , / 2 /<t 1+ l ' V 


+ (v TE ) i+l / 2 / ’ V z ~ &Zp 2 ^ rv p^ 2 ,r “ * 2 ,z TE ^i+l /2 


+ Az 0 o c frK/g (v 4 > 0 ) (v V )]. )/0 

2 <j> p 2 ,r ,r p z,r i+ 1/2 


6.15 


For $ 2 : 


A t = max {0, (DF$ - 0.5 CONj) /2 , - C 0 N$ i+1/2 ] , 


6.16 


B i = max [0, (DF$ + 0.5 CON® ) 1—1 / 2 ^ * 


6.17 


C l-»2 1S V2 , l ’ 


6.18 


D i -»2 ' S Vl +A 1-1 + B l + 1 + D^<AS,I 


6.19 


where 


CO'S - »2 <r Vi.l/2 ( “ 2 - IE) i 


6.20 


D». -p 2 4 Z (to /o ) i+1/2 U/(r i+1 - r x ) 


TV^l (1 ■S I “ +1 /» 2 I 1+1 > . 


6.21 


Note that the A i and formulae (6.4, 6.5, 6.10, 6.11, 6.16 and 6.17) 

are hybrid in nature to account for high lateral convection (see Spalding, 
~2 

1978). v in Equations 6.9 and 6.15 is calculated from Equation 4-10. 
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The equation set generated by Equation 6.3 for the different nodes at any 
station is solved by using the Tri-Diagonal Matrix Algorithm (TDMA) to obtain 
the dependent variable at the different nodes (see Roache, 1972). 

At the beginning of the calculations, ^ is unknown. Thus, the C0MP2P 
subroutine uses a guess-and-correct procedure for 4^ or The procedure of 

correcting 4^ or 4^ takes place in the subroutine ADJ2P. There are three 
procedures for correcting 4^ : direct substitution, computed under relaxation, 

and use of pressure corrections. In the direct substitution procedure the 
computed values of 4^ are used as the guessed ones for the start of the new 
solution loop. The disadvantage of this procedure is the possible non- 
convergence, when the changes are large. The advantages of this method are 
its simplicity and the consequent economy. This method can be expected to 

_3 

work satisfactorily when 4^ values are small ( 4„ <10 ). 

The computed under-relaxation method starts the next iteration loop 
with given by- 

*1 " *1 + 5 < * 1 '“ * 1 >* 6.22 

where the * denotes old values and X, is an under-relaxation factor, 
conveniently taken as: C a 4^ | The under-relaxation factor is expected 

to be slight for dilute suspension, i.e., this method reduces to the direct 
substitution method. 

The pressure corrections method devises and solves a pressure correction 
equation driven by errors: (1 - ^ ~ ^). This method is suitable for' large 
values of 9 , which are outside the scope of the present study. 
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6.3 The Solution Procedure 


The steps to obtain the solution at a given axial location are: 

1. Guess the downstream <J> ^ distribution (from the upstream values). 

s s 

2. Solve for U downstream: obtain r* and U * . 

z r 

3. Solve for K and e to obtain the eddy diffusivities. 

4. Calculate the local dimensionless quantities, Re^, Sc^, and Sh^. 

Then obtain the mass transfer rate; m , hence the sink terms in the mass 
conservation equation for each group or the total source tem in the 
continuity equation for the carrier phase. 

5. Obtain the downstream diameter distribution from the upstream values and 
the local evaporated mass. 

6. Calculate the size range for each group from knowing the largest and 
smallest droplet diameters, and the number of sizes to be considered (it 
could be different from the upstream value). 

7. Label the droplets according to their local diameters and the size 
ranges for each group. 

8. Obtain p(r) from the gas-phase lateral momentum equation. 

Ic Ic Ic If * c 

9. Solve for downstream (v /V ) , V and a and get 4,. <f 

r z ’ z z 1 

* 

10. Compare the new $ ^ with the guessed $ ^ . 

11. Make corrections and repeat steps 1-7 until the solution converges 
before marching to the next station. 

In general, two iterations are needed at each station to achieve convergence. 

6.4 The Boundary Conditions 

The parabolic flows considered in this report require the prescription of 
three boundary conditions for each dependent variable. At the predictions 
starting plane the profiles of all the dependent variables must be specified 
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from the experimental data. At the axis of symmetry (r = 0) all the radial 
gradients are set to zero, in. addition to the vanishing radial velocity. The 
jet boundaries are determined via the adjustment of the entrainment rate and 
the specification of the radial gradient of U z to fixed small value. Just 
outside these boundaries the values of the other dependent variables will be 
those corresponding to the ambient conditions. For example, all V 's 
and are ec l ua l to zero there.. 

In the next section the results presented are obtained- using 40 lateral 
nodes to span the flow domain between the centerline of the jet and its outer 
edge. Grid-dependence tests were conducted with 30, 40, and 50 lateral nodes 
and different axial step sizes and it was concluded that the 40 node grid 
results are virtually grid-independent. 
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7.0 RESULTS 


First, a turbulent round gaseous jet laden with monosize solid particles 
is considered. This flow allows the study of the interaction between the two 
phases and the particles dispersion due to the turbulence without the 
complexity of mass transfer. The predictions are compared with the data of 
Modarress et al . (1984) in section 7.1. Since no experimental data exists for 
an evaporating spray in the developing region (z/D < 20), the model is 
considered to predict an idealized flow of a turbulent round jet laden with 
multisize evaporating methanol droplets in section 7.2. Two more cases are 
considered in sections 7.3 and 7.4. Both of these two cases are for a Freon- 
11 spray issuing from an air atomizing nozzle where experimental data are 
available. The first flow is that of Shearer et al . (1979) where the data are 
available at distances equal to or greater than 170 nozzle diameters. The 
second flow is that of Solomon et al . (1984) where the data are available at 
distances equal to or greater than 50 nozzle diameters. In both cases the 
predictions are compared with the measurements. Table C-l (Appendix C) 
summarizes the different cases considered in this study. 

7.1 The Flow of Modarress et al . (1984) 

Modarress et al. (1984) reported much needed experimental data to help 
understand the behavior of two-phase turbulent jets and validate their 
theoretical models. They investigated the effects of 50 ym and 200 ym glass 
particles on the mean air velocity and the turbulent stresses at two different 
mass loading ratios, 0.32 and 0.85. Figure 1-2 shows a sketch of the two- 
phase turbulent jet considered by Modarress et al . (1984). Air carrying 
uniform-size glass particles issues vertically downwards from a cylindrical 
pipe of diameter D, 0.2 m. The jet is enclosed in a cylindrical container 
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with a diameter equal to 30 D to avoid ambient disturbances. An air stream of 
low velocity surrounds the nozzle and extends to the container wall to provide 
the required entrained mass by the jet, thus preventing the occurrence of 
internal circulation in the measurements region. Table C-2 (Appendix C) lists 
the experimental conditions at 0.1 D downstream of the pipe exit. These 
values represent the initial conditions for the dependent variables required 
in the numerical calculations. 

In what follows the predicted distributions are compared with the 
measured distributions of the mean velocities, volume fractions of the two 
phases, turbulence intensity and shear stress of the gaseous phase and the jet 
spreading rate. Figures 7-1 and 7-2 show the effects of the particles; mass 
loading ratio (X Q = 0.32 versus 0.85) on the mean velocities for 50 p 
particles (Case 1 and 2). Figures 7-2 and 7-3 show the effects of the 
particles; diameter (50 P versus 200 p) at almost the same loading ratio (0.8) 
on the mean velocities of the two phases. First the main effects of the 
particles on the carrier phase velocity are discussed, then the behavior of 
the particles' velocity and volume fraction at the different mass loading 
ratios are discussed. 

Figures 7-1 to 7-3 show the radial profiles of the mean axial velocities 
of the two phases at z/D = 20, normalized by the corresponding mean centerline 
velocity of the single-phase jet, U z c . s * The flow conditions are those of 
Cases 1, 2 and 3 in Table C-2 (Appendix C). Also shown is the mean velocity 
profile of the turbulent single-phase jet having the same Reynolds number 
(14100) at the pipe exit. 

It can be seen from these figures that the mean velocity of the carrier 
phase is highly affected by the presence of the particles in the inner region, 
especially at the jet centerline (30% higher than that of the single phase for 
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FIGURE 7-1 RADIAL DISTRIBUTION OF THE MEAN 
AXIAL VELOCITIES AT z/D ■ 20 (CASE I ) 
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FIGURE 7-2 RADIAL DISTRIBUTION OF THE MEAN 
AXIAL VELOCITIES AT i/D » 20 (CASE 2) 


118 



FIGURE 7-3 RADIAL DISTRIBUTION OF THE MEAN 
AXIAL VELOCITIES AT z/D-20 (CASE 3) 
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Case 1, 75% for Case 2, and 25% for Case 3). This behavior can be explained 
by the fact that particles are confined to the inner region of the jet. Due 
to this confinement and the high inertial forces of the particles, their 
centerline velocity decays with the downstream distance at a slower rate than 
that of the fluid (see Figure 7-8), and thus they become a source of momentum 
to the fluid. Also due to the confinement, the number density of the 
particles which is a strong parameter in the momentum transfer between the 
particles and the carrier phase at any cross section is maximum at the jet 
centerline. The “conf inement of the particles is evident in Figures 7-2 and 7- 
3 where the concentration of the solid particles vanishes at a radial distance 
of r/z = 0.06, while the fluid spreads to at least three times this 
distance. This confinement can be explained by the fact that heavy particles 
do not respond well to fluid turbulence fluctuations ( << v t ), thus the main 
force that accelerates the particles in the radial direction is the viscous 
drag. This drag force is proportional to (U r -V r ), and since U r is negative in 
the outer region of the jet and V r < U r ), the resulting force will be directed 
inwards thereby limiting the radial spread of the particles. 

The influence of the loading ratio of the dispersed phase on the mean" 
velocities at z/D = 20 for 50 p particles is displayed in Figures 7-1 and 7- 
2. The inlet conditions for the two cases (1 and 2) are identical except for 
the loading ratio. By increasing the loading ratio from 0.32 (Case 1) to 0.85 
(Case 2) the carrier phase velocity at the centerline increases from 30% to 
75% relative to the corresponding velocity of the single-phase jet. It can be 
seen also from those two figures that the ratio between the centerline 
velocity of the dispersed phase to that of the single phase is 1.5 for Case 1 
and 2 for Case 2. This can be explained by the fact that the initial momentum 
of the dispersed phase is proportional to the mass flow rate of that phase 
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since they have the same velocity distribution. Therefore, the initial 
momentum of the dispersed phase for Case 1 is about 2.7 times that for Case 
2. This enhances the momentum transfer from the particles to the air, thereby 
enhancing the increase in the axial velocity of the latter compared with that 
of the single phase. Now, since the axial air velocity at any point on the 
jet axis at the higher loading case is greater than the corresponding value at 
the lower loading case, the momentum drain from the particles is expected to 
be inversely proportional to the particles' loading ratio. This explains the 
higher centerline velocity of the particles at the higher mass loading ratio. 

The influence of the particles' diameter of the dispersed phase on the 
mean velocities at z/D = 20 is displayed in Figures 7-2 and 7-3 (Cases 2 and 
3). The main difference between the two cases is the particle diameter, so 
any quantitative change in the mean velocity profiles is attributed to two 
factors: 1) the interphase surface area or the momentum exchange coefficient, 
and hence the source terms of the momentum equations, and the K and e 
equations; 2) the particles' response to the turbulent fluctuations, thus the 
additional turbulence dissipation caused by the fluctuating particle slip 
velocity and its correlation with the fluid velocity fluctuation. The surface 
area in Case 2 is four times that in Case 3, since, for nearly the same 
loading ratio, the number of the 50 y particles is 64 times that of the 
200 y particles. This increase in the number of particles or the interphase 
area results in augmenting the momentum source of the carrier fluid 
consequently reducing the rate of decay of its centerline velocity (see Figure 
7-8). 

Figure 7-4 shows the reduction in the shear stress due to the existence 
of the particles with the air in the same control volume for 50 y and 
200 y particles. Due to the reduction in the turbulence kinetic energy 
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FIGURE 7-4 RADIAL DISTRIBUTION OF THE 
SHEAR STRESS AT z/D-20 ISOLID PARTICLES) 
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compared with that of the single phase, which is associated with additional 
turbulence dissipation, the carrier phase momentum diffusivity is reduced. 
Figure 7-4 shows that the reduction in the shear stress for 50 u particles is 
greater than that of 200 p particles. This reduction in the turbulent 
diffusion coefficient reduces the rate of decay of the centerline air mean 
velocity (see Figure 7-8). 

It is also clear from Figures 7-1, 7-2, and 7-3 that the single-phase jet 
is wider than the particle-laden jet; this will be discussed later in this 
section. Figures 7-1, 7-2, and 7-3 display in general good agreement between 
the measured and predicted velocity and concentration profiles. 

In order to distinguish between the dispersed phase effects on the mean 
motion (inertia and drag) and on turbulence (diffusion), the mean velocity 
profiles obtained by solving the governing equations for the mean motions 
(Equations 3.26 to 3.33 together with the single-phase K and e equations 
(i.e., Equations 3.34 and 3.35 without the additional production and 
dissipation terms due to the dispersed phase) are shown in Figure 7-1 and 7- 
2. The resulting increase in the fluid centerline velocity, as compared with 
that of the single-phase jet, is only half that measured and predicted by the 
new K-e model. Stated differently, the modulation of the fluid mean velocity 
profile by the dispersed phase is not only due to the particles' Inertia and 
drag but equally due to the additional turbulence dissipation. This in turn 
reduces the fluid momentum diffusivity resulting in a peaked velocity profile 
near the jet centerline. The additional turbulence dissipation is caused 
mainly by the fluctuating particle slip velocity and its correlation with the 
fluid velocity fluctuation that appeared in Equations 3.34 and 3.35. The 
consequent reduction in the fluid turbulence intensity and shear stress is 
displayed in Figures 7-4 and 7-5 where the agreement between the measurement 
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and prediction is good. 

Figure 7-6 compares the concentration distribution of 50 y particles 

(Case 2) with that of 200 y (Case 3) particles. Since the mass eddy 

2 -2 

diffusivity is inversely proportional to (U-V) /v (Equation 5.30), it is 
consequently higher for the 200 y than that of the 50 y particles, so one 
would expect that the particles of 50 y will diffuse in the radial direction 
more than the of 200 y. This is evident in Figure 7-6 where the agreement 
between the measurement and prediction is good. 

Figure 7-7 shows the effect of the dispersed phase on the spreading rate 
of the jet by comparing the different “ z distributions of the three 

cases, where Yj /2 *- s th e radius at which the fluid mean axial velocity is half 
that at the centerline. While for a turbulent single-phase jet the value of 
the slope (dYj^/^ 2 ) is constant (= 0.08), the value for a two-phase jet is a 
function of the dispersed phase properties such as particle diameter and 
density and loading ratio. This dependence is displayed in Figure 7-7. For 
Case 3 (d = 200 y, X q = 0.8) the predicted slope value is 0.053, for Case 2 

(d = 50 y, X q = 0.85) it is 0.046, and for Case 1 (d = 50 y, X Q = 0.32) it is 

0.064. Cases 3 and 2 ' have nearly the same loading ratio but the particle 

diameter in the latter is one quarter that of the former, the result being a 

reduction of the spreading rate by more than 13%. 

Figure 7-7 also shows the discrepancy that results in predicting the 
spreading rate if the single-phase K-e model is used instead of the proposed 
model. The former predicts for Case 1 a slope of 0.072 while the latter 
agrees with the experimental value of 0.064. As explained earlier this is due 
to the fact that the additional dissipation of turbulence energy as a result 
of the dispersed phase is accounted for in the proposed model. 

Figure 7-8 shows the decay of the mean centerline velocities of the two 
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FIGURE 7-6 RADIAL DISTRIBUTION OF THE 
PARTICLES VOLUME FRACTION AT x/D»20 
(SOLID PARTICLES) 
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phases for Cases 1 and 3 compared with the single-phase values. Here U z Q is 
the carrier-phase centerline velocity at the pipe exit. It can be seen from 
Figure 7-8 that the two phases reach a local equilibrium situation, equal 
velocities, at about 10 pipe diameters and after that the relative velocity 
between the two phases along the jet centerline increases by increasing the 
downstream distance from the pipe exit. This behavior was pr ev i° us ly analyzed 
in the discussion of Figures 7-1 through 7-3. 

7.2 The Methanol Spray 

The flow considered in the present study is identical to the flow of 
Modarress et al. (1984) except that the solid spheres are replaced by methanol 
droplets of a given size distribution at the exit of the pipe (Figure 1-2). 

The goal here is to mimic the flow of an idealized spray that has well-defined 
initial conditions. In the present study the good agreement between 

prediction and experimental data in the cases of a round gaseous jet laden 

with solid particles allows the use the of the latter while adding the 

complexity of mass transfer and the resulting size changes in the same jet. 

A turbulent round jet laden with multisize evaporating liquid droplets is 
considered in this section. Atmospheric air carrying methanol liquid 
droplets of diameters 100, 80, 60, 40 and 20 tfn issues vertically downwards 
from a cylindrical pipe of diameter D (= 0.02 m) . The initial mean velocity 
and the turbulence intensity distributions are assumed to be those of the 
fully developed pipe flow as in the work of Modarress et al. (1984). The 
ratio between the velocity of the dispersed phase to that of the carrier phase 
at the centerline is equal to 0.7. The carrier fluid Reynolds number is equal 
to 30,000. The temperature of methanol droplets is assumed to be uniform at 
the steady state saturation conditions. The initial mass flow rates of the 
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different size groups are assumed to be equal and have a plug profile for 
volume fractions. Three different mass loading ratios of 0.1, 0.25, and 0.5 
(Case 6) are considered. 

In what follows, the predicted mean velocities, volume fractions of each 
phase, turbulence intensity, and shear stress of the carrier phase under the 
three mass loading ratios are presented. 

The normalized radial profiles of the mean quantities of the different 

phases at 20 pipe diameters from the exit plane at X Q = 0.5 are shown in 

Figure 7-9. The mean velocities of the carrier phase and those of the five 

groups (k = 1,2,..., 5) of droplets are normalized by the centerline velocity 

of the single phase jet, U z CiS .« Here k = 1 refers to the group that has the 

largest diameters, and k = 5 the smallest ones. It can be seen from this 

figure, as one expected, that the difference between the velocity of the 

carrier phase and that of the largest diameter group is greater than that of 

any other group. This is attributed to the balance between the inertia of the 

droplet and the momentum exchange’ force. The inertia force is proportional to 

(d^)^ whereas the momentum exchange force is proportional to the droplet 

diameter with an exponent ranging from 1 to 1.7 (for a Reynolds number less 

than 100). If all the turbulent correlations in Equation 3.31 due to their 

small values compared with the mean momentum exchange term are now dropped, 

lc 

the equation becomes independent of $ . If the droplet size is then Increased, 
the inertia becomes much greater than the momentum exchange force, and as a 
result the relative velocity between the droplets and the carrier phase (U z - 
V z k ) increases. The volume fraction profile of each group normalized by the 
centerline value of the first group is shown also in Figure 7-9. Since the 
reduction rate of the droplet diameter due to the evaporation process is 
inversely proportional to the square of the diameter, the smaller the droplet 
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FIGURE 7-9 RADIAL DISTRIBUTION OF THE NORMALIZED 
MEAN VELOCITIES AND VOLUME FRACTIONS AT z/D= 20 
AND AT X e ■= 0.5 (CASE 4 ) 
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diameter is, the more reduction in the volume fraction. Figure 7-9 also shows 

that the smaller the mean droplet diameter is, the less peaked the volume 

fraction profile of its group. This is attributed to the turbulent diffusion 
k 

coefficient ( v ) of the droplet which decreases with the increase of the 
P 

relative velocity or the droplet diameter (Equation 5.30). It can also be 
seen from Figure 7-9 that the mean velocity of the carrier phase is affected 
by the presence of the dispersed phase especially in the inner region. 
Elghobashi et al. (1983) discussed in detail how the entrainment and the 
negative radial velocity of the carrier phase in the jet outer region 
influence the volume fraction distribution of the dispersed phase. They 
showed that the entrainment flow creates an inward force exerted on the 
droplets towards the jet centerline. This force combined with the small 
turbulent diffusivity of the droplets, compared with that of the carrier 
phase, renders the volume fraction profile of the dispersed phase 
significantly narrower than the velocity profile of the carrier phase. Since 
the momentum exchange between the two phases is a linear function of the 
droplets volume fraction, it could be expected, that the momentum transfer to 
the carrier phase is maximum at the jet centerline. At the same time the 
reduction in the turbulence kinetic energy of the carrier phase and the 
increase of the dissipation rate of that energy due to the presence of the 
dispersed phase in the same control volume lead to a less turbulent diffusion 
coefficient for the carrier phase and hence a less radial diffusion of that 
phase compared with the single phase. These two factors make the velocity of 
the carrier phase at the jet centerline much greater than that of the single- 
phase jet (30% higher) and less than its corresponding value in the jet outer 
region. 

The Influence of the loading ratio of the dispersed phase on the carrier 
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fluid turbulence intensity and shear stress is displayed in Figures 7-10 and 
7-11. The reduction in the turbulence energy or the increase in the 
dissipation rate of that energy is caused by the fluctuating relative velocity 
between the droplets and the carrier phase and the turbulent correlation 
between this velocity and other fluctuating quantities, volume fractions and 
carrier fluid velocity. It can be stated that the reduction in the turbulence 
intensity and the shear stress is proportional to the mass loading ratio but 
not linearly. 

The concentration of the evaporated material in the carrier phase is 
shown in Figure 7—12 at two different axial locations (z/D = 10 & 30) and at 
X Q = 0.5. Due to the continuous air entrainment by the jet and the turbulent 
diffusion of the vapor, the concentration of the evaporating material in the 
carrier fluid at z/D = 30 is less than the corresponding values at z/D = 10 at 
the same distance from the jet axis, although the total evaporated mass 
increases with downstream distance. This is also true even at the jet 
centerline as will be seen in the discussion of Figure 7-15. 

It can be seen from Figure 7-12 that C is minimum in the jet outer region 
and maximum at the jet centerline. Since C^, according to the assumption of 
this study and the droplets' material, has a constant value of 0.12, the 
transfer number (Equation 3.42) is maximum in the outer region of the jet. 
Therefore, the diminution rate of the droplet diameter is greater in the outer 
than in the inner region. 

Figure 7-13 shows the centerline decay of the mean axial velocities of 
the different groups and the carrier phase compared with the single phase 
values for Xq = 0.5. Here U z Q is the carrier-phase centerline velocity at 
the pipe exit. It can be seen that the relative velocity between the droplets 
and the carrier phase or the disequilibrium of the flow along the jet 
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FIGURE 7-10 RADIAL VARIATION OF THE TURBULENCE 
INTENSITY WITH X 0 AT z/D=20 ( CASE 4 ) 
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FIGURE 7-11 RADIAL VARIATION OF THE SHEAR 
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FIGURE 7-12 RADIAL DISTRIBUTION OF THE 
VAPOR CONCENTRATION (CASE 4 ) 
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centerline, increases with increasing the droplet's diameter. It is worth 
noting that the carrier-phase centerline velocity is about 30% higher than the 
corresponding value of the single phase in the range 7 < z/D < 30 as 
previously discussed. 

Figure 7-14 exhibits the centerline decay of the volume fraction and mean 
droplet diameter based on the total surface area of the droplets for the five 
groups. The mean diameter is a quantity that is not used in any calculations 
but facilitates the display and discussion of the results. In the present 
work, it was possible to calculate the local diameter distribution within each 
group, thus from the maximum and minimum diameters at any station and the 
number of sizes to be solved, the diameter range for each group can be fixed 
(e.g., at z/D = 10, group k = 1 contains droplets ranging from 95 to 78 
microns). It can be seen from Figures 7-9 and 7-14 that the smaller the 
droplet diameter is, the higher the evaporation rate, hence the rapid decay of 
the volume fraction and the mean diameter. 

Figure 7-15 shows the axial distribution of the total volume fraction of 

the droplets and the centerline concentration of the methanol vapor in the 

carrier phase (C ") for different mass loading ratios. Here <t_ /$» is the 

* Z j C Z j o 

total volume fraction of the dispersed phase at the centerline divided by that 
value at the pipe exit. The concentration of the evaporated material in the 
carrier phase first increases until z/D = 10 then monotonically decreases due 
to the continuous air entrainment by the jet and turbulent diffusion of the 
vapor. 

The variation of the maximum turbulence intensity and maximum shear 
stress of the carrier phase with the axial distance is displayed in Figures 7- 
16 and 7-17 for the different mass loading ratios. It can be seen that the 
reduction in the turbulence quantities is proportional to the mass loading 
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FRACTION AND VAPOR CONCENTRATION WITH X, ( CASE 4 ) 
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FIGURE 7-16 AXIAL VARIATION OF THE MAXIMUM 
TURBULENCE INTENSITY WITH X fl ( CASE 4 ) 
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FIGURE 7-17 AXIAL VARIATION OF THE MAXIMUM 
SHEAR STRESS WITH X„ ( CASE 4 ) 
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ratio but again not linearly. These two figures also show that farther 
downstream from the pipe exit, the turbulence quantities are approaching their 
values for a single-phase jet due to the continuous diminution of the 
droplets' volume fraction. 

The rate of evaporation is a function of both the transfer number and 
droplet Reynolds number, which are maximum in the outer region and minimum at 
the centerline. So the rate of evaporation is maximum in the outer region or 
the minimum droplet diameter. This explains the radial distribution of the 
droplet diameter at the various sections as shown in Figure 7-18. Also 
displayed is the monotonic reduction in droplet diameters with distance 
downstream for the five groups. 

Figure 7-19 shows the effect of the evaporating spray on the spreading 
rate of the jet by comparing the different ~ 2 distribution, where Yj^ 

is the radius at which the carrier-fluid mean axial velocity is half its value 
at the centerline. While for a turbulent-single phase jet the value of the 
slope (dYj^/d 2 ) * s constant (= 0.08), that for a two-phase jet is a function 
of the dispersed phase properties such as droplet diameter, density and mass 
loading ratio. This dependence was discussed in the work of Mostafa and 
Elghobashi (1983). In the developing region, the spreading rate of the spray 
case is much less than that of the single phase. As vaporization proceeds the 
effects of the droplets on the carrier fluid diminish allowing the fluid 
behavior to approach that of a single-phase jet. 

7.3 The Flow of Shearer et al. (1979) 

Shearer et al. (1979) measured the carrier phase properties in a 
turbulent two-phase round jet using a laser doppler anemometer, the droplet 
size distribution and the liquid mass flux using the inertial impaction 
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method. The Freon-11 spray was generated by an air-atomizing nozzle of 1.194 
mm outer diameter (D). The ratio of the mass flow rate of Freon-11 spray at 
the nozzle exit to that of the air (X Q )~ is equal to 6.88 and the initial 
average velocity, U z 0 = 74.45 m/s. They also measured the mean mixture 
fraction by isokinetical ly sampling the flow at the gas velocity. Shearer et 
al. (1979) measured the radial profiles of the mean and rms velocity, and the 
Reynolds stress at three stations (z/D = 170, 340, and 510) for both 
isothermal single- phase and vaporizing spray jet flows. For computational 
purposes, the profiles of turbulence dissipation rate ( e) at z/D = 170 are 
obtained from the shear stress measurements and the axial velocity gradient at 
the same station (z/D = 170). Also, the velocity distribution of the droplets 
(one group with an average diameter = 27 pm) is assumed to be the same as that 
of the carrier phase. This assumption will be discussed at the end of the 
next section. From the measurements of the droplets' mass flux and velocity 
distribution, the volume fraction (^) is obtained. The profile of the freon 
vapor concentration in the carrier phase (C) is obtained from the mixture 
fraction measurements and the state relations given by Shearer (1979). Table 
C-3 (Appendix C) summarizes all the starting profiles needed for the 
computation for both the single-phase jet and the evaporating spray cases. 

Temperature measurements of the carrier phase (with a bare wire 
chromelalumel thermocouple) showed only 5°C difference either in the radial or 
the axial directions (between z/D = 170 & 510). On the other hand, Shearer's 
analysis (1979) showed that the droplet's temperature at z/D = 170 is equal to 
the Freon's saturation temperature (240. 3°K). In the present calculations it 
was assumed that the temperature of the carrier phase is equal to the 
surrounding air temperature (296°K) and the droplet's surface temperature is 
equal to the saturation one (240. 3°K). At these conditions, the density of 
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the liquid Freon-11 is equal to 1518 Kg/ra^ and the vapor concentration at the 
droplet's surface (C^ in Equation 3.42) is equal to 0.292. In what follows we 
compare the predicted with the measured distributions of the mean velocity and 
shear stress of the carrier phase at z/D = 340 and 510. 

Figure 7-20 shows the measured and predicted centerline decay of the mean 
axial velocity of the carrier phase compared with the single-phase values. 

Due to the fact that the inertia of the droplets is much greater than that of 
the carrier phase (p^/ = 1500) , the centerline velocities of the droplets 
in the region z/D < 170 are greater than those of the gas. As a result, the 
centerline velocity of the carrier phase would be expected to be greater than 
that of the single phase. This is due to 1) the continuous momentum transfer 
from the droplets to the gas since V z c is greater than U z c in the region 
close to the nozzle (z/D < 170) and 2) the reduction of the turbulence 
intensity (and hence turbulent diffusion) in the spray case compared with that 
of the single-phase jet (as will be seen later in Figures 7-23 and 7-24). 

Figures 7-21 and 7-22 show the normalized radial profiles of the mean 
axial velocities at 340 and 510 nozzle diameters from the exit plane for both 
the single-phase jet and the evaporating spray cases. It can be seen from 
these figures that the jet width in the spray case is narrower than the 
single-phase one. This result can be attributed to the increase of the 
centerline velocity of the carrier phase compared with its corresponding value 
in the single-phase jet. The experimental data show that with increasing the 
distance downstream from the nozzle exit, the jet behavior approaches that of 
the single phase (Figure 7-22). The effect of the droplets on the radial 
shear stress distribution is displayed in Figures 7-23 and 7-24. It should be 
noted that the starting values of the turbulence quantities and mean velocity 
distribution of the vaporizing spray case differ (less shear stress) from 
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FIGURE7-20 AXIAL VARIATION OF THE MEAN CENTERLINE VELOCITY (CASE 5 ) 
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FIGURE 7-21 RADIAL VARIATION OF THE MEAN AXIAL 
VELOCITY AT z/D = 340 (CASE 5 ) 
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FIGURE 7-22 RADIAL VARIATION OF THE MEAN AXIAL 
VELOCITY AT z/D = 510 (CASE 5 ) 
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FIGURE 7-23 RADIAL VARIATION OF THE SHEAR STRESS 
AT i/D = 340 (CASE 5 ) 
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FIGURE 7-24 RADIAL VARIATION OF THE SHEAR 8TRESS 
AT z/D = 510 (CASE 5 ) 
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those of the single phase. This may have some effects on the profiles 
downstream (z/D = 340). In general, there is a reduction in the shear stress 
or an increase in the dissipation rate of the turbulence kinetic energy due to 
the presence of the liquid droplets in the same control volume with the 
carrier phase. As vaporization proceeds the effects of the droplets on the 
turbulence quantities diminish allowing the fluid behavior to approach that of 
a single-phase jet (Figure 7-24). 

In the present case it was assumed that the velocities of the droplets, 
are equal to those of the gas. To study the effect of this assumption on the 
results, the droplets' velocity was increased by 20%; the effect on the 
carrier phase profiles was negligible. This result can be attributed to two 
factors: 1) the droplets' diameter, at the starting station, is equal to 
27 Pm or less, so the reduction rate of the mean slip velocity between the 
droplets and gas is considerable due to the vaporization; 2) since the 
droplets' mass fraction was measured, an increase in the velocity necessitates 
a reduction in the volume fraction. Thus, the effects of increased velocity 
are counterbalanced by those of decreased volume fraction. 

It is important to note that the effects of density fluctuation in the 
calculation were neglected here. This assumption can be justified in the 
present study since the mean density gradient is very small compared with the 
velocity gradient. This is due to the negligible evaporated mass compared 
with the entrained air, so the properties of the carrier phase are almost 
those of the standard air. 

7.4 The Flow of Solomon et al . (1984) 

Solomon et al . (1984) presented some comprehensive measurements of the 
detailed structure of a two-phase turbulent round jet. Experiments considered 
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the same test rig of Shearer et al . (1979) to perform some detailed 
measurements of the droplets' properties. Two mass loading ratios (X Q = 7.71 
and 15.78) were considered. Solomon et al . measured the carrier phase 
properties using a laser doppler anemometer, the droplet size and velocity 
using the shadow photograph technique, the liquid mass flux using inertial 
impaction method, and the mean mixture fraction by isokinetically sampling the 
flow. The radial distributions of these quantities were reported at four 
stations (x/D = 50, 100, 250, and 500) for the two mass loading cases. They 
classified the droplets into finite-size groups and measured the velocities 
and the number density distribution of each group. For computational 
purposes, the profiles of the turbulence dissipation rate ( e) , the volume 
fraction of each droplets group ( $ ) , and the freon vapor concentration in the 
carrier phase ( C) are obtained from the different measured quantities at z/D = 
50. e is obtained from the distributions of the turbulent shear stress, the 
mean axial velocity gradient, and the turbulence kinetic energy at the same 
station. Seven groups (17.5, 22.5, 27.5, 32.5, 42.5, and 52.5 pm) are 
considered for X Q = 7.71 (Case 6) and ten groups (15, 25, 35, 45, 55, 65, 75, 

If 

85, 95, and 100 Pm) are considered for X Q = 15.78 (Case 7). 4> is obtained 

from the distributions of the liquid mass flux, and the mean velocity of the 
different droplets groups and their relative number density at z/D = 50. C is 
obtained from the mixture fraction measurements and the state relations given 
by Solomon et al. (1984). Table C-4 (Appendix C) summarized all the starting 
radial profiles of the main dependent variables at z/D = 50. This information 
is essential for accurately predicting the present flow to validate the 
turbulence model put forth in the present study. 

Temperature measurements of the carrier phase (with a bare wire 
chromelalumel thermocouple) showed a maximum temperature difference of only 
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20°C either in the radial or in the axial directions. The analysis of Solomon 
et al. showed that the droplet's temperature reaches the Freon's saturation 
temperature at z/D = 50. It was assumed in the present calculations that the 
temperature of the carrier phase is equal to that of the surrounding air 
(300°K) while the droplet surface temperature is equal to that of the 
saturation conditions (240. 3°K). 

A comparison of the predicted with the measured distributions of the mean 
velocity, turbulence intensity and shear stress of the carrier phase, and the 
mean velocity of the droplets of the different-size groups follows. 

. Figures 7-25 and 7-26 show the measured and predicted centerline velocity 
distributions of the carrier phase and those of the different droplet groups 
for cases 6 and 7. Here k = 1 refers to the group that has the largest 
diameters, and k = 7 or 10 the smallest ones. The mean velocities are 
normalized by the average velocity at the nozzle exit (U Q = 64.5 m/s for case 
6 and 29.64 m/s for case 7). It can be seen from the figures that the 
relative velocity between the carrier phase and the group of largest diameters 
is greater than that of any other group. This behavior is already explained 
in the analysis of Figure 7-9 (section 7.2). Figures 7-25 and 7-26 also show 
the continuous reduction in the relative velocity between the carrier phase 
and the group of smallest diameters as the distance measured from the nozzle 
exit plane increases. This can be attributed to the fact that the smaller the 
droplet diameter, the higher the reduction rate in the droplet diameter 
itself. Thus, by increasing the downstream distance, the smallest diameters 
group satisfies the local equilibrium conditions where the velocity of the 
droplet becomes equal to that of the carrier phase. It can be seen also from 
Figures 7-25 and 7-26 that the relative velocity between the carrier phase and 
the group of largest diameters increases with an increase in the downstream 
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FIGURE 7-26 THE DISTRIBUTION OF THE JET CENTER LINE MEAN VELOCITIES 
(CASE 7 ) 



distance. This behavior is already explained in the analysis of Figure 7-13 
(section 7.2). Figures 7-25 and 7-26 display in general good agreement 
between the measured and predicted mean centerline velocities. 

The influence of the loading ratio of the dispersed phase on the 
centerline mean velocity distributions of the carrier phase is illustrated in 
Figure 7-27. In this figure the increase in the mean centerline velocity of 
the carrier phase compared with the corresponding value of the single-phase 
jet is proportional to the mass loading ratio (but not linearly). This 
proportionality is analyzed in detail in the discussion of Figures 7-1 and 7-2 
in section 7.1. Figures 7-28 to 7-30 show the normalized radial profiles of 
the mean axial velocities of the carrier phase at 100, 250, and 500 nozzle 
diameters from the nozzle exit plane for both the two loading ratios (Cases 6 
and 7). It can be seen from these two figures that the jet width decreases 
with the increase of the mass loading ratio. 

Figure 7-28 shows a maximum discrepancy of 30% between the predicted and 
measured velocities although the agreement is very good in Figures 7-29 and 7- 
30. Probably the measured quantities are overestimated at this station since 
Solomon et al. reported the same discrepancy between the measurements arid 
their predictions at the same station, using the Lagrangian frame of work. 

The influence of the dispersed phase on the carrier fluid turbulence 
kinetic energy and shear stress is displayed in Figures 7-31 to 7-36. It can 
be stated that the reduction in the turbulence energy and the shear stress is 
proportional to the mass loading ratio but not linearly. These figures also 
show that farther downstream from the nozzle exit (z/D = 500), the turbulence 
quantities are approaching their corresponding values for a single-phase jet 
(based on the experimental data of Shearer et al., 1979). 

To understand the nature of the turbulent interaction between the carrier 
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FIGURE 7-28 THE RADIAL DISTRIBUTION OF THE 
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FIGURE 7-29 THE RADIAL DISTRIBUTION OF THE 
MEAN GAS VELOCITIES AT z/D«250 
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FIGURE 7-30 THE RADIAL DISTRIBUTION OF THE 

MEAN GAS VELOCITIES AT z/D-500 
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FIGURE 7-33 THE RADIAL DISTRIBUTION OF THE 
KINETIC ENERGY OF TURBULENCE 
AT i/O* 500 
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fluid and the vaporizing droplets, the main features of this type of flow are 
summarized as follows: 

1. The toal mass of the dispersed phase continuously decreases and so does 
the volume fraction. Due to the reduction of the volume fraction, the 
momentum exchange terras (mean and/or fluctuating) are reduced. 

2. The velocity of the evaporating material as it leaves the droplet 
surface is different from that of the carrier fluid. Thus, there is an 
additional momentum transfer that depends on the evaporation rate and the 
relative velocity. 

3. The momentum exchange coefficient is inversely proportional to the 
droplet diameter with an exponent ranging from 2 to 1.3 (for a Reynolds 
number less than 100). Hence, as the diameter is reduced the momentum 
exchange coefficient increases. 

4. The vaporization reduces the droplets* diameter and thus the total 
relative mean velocity (U - V) and the higher the turbulent diffusivity of 
the dispersed phase is. 

Figures 7-31 to 7-36 display in general good agreement between the 
measured and predicted turbulence kinetic energy and shear stress of the 
carrier phase. 

The predictions of the axial distribution of the Sauter mean diameter at 
the jet centerline compared with the experimental data is displayed in Figure 
7-37. This diameter is given by 


i (d ) 3 n /V 

SMD = - 7.1 

i ± (<V n /v 


where n^ is the number of the droplets of diameter d^. It is clear that there 
is also" good agreement between the predictions and the data for the averaged 
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diameter. 


It is important to note that the present study neglected the effects of 
density fluctuations in these two cases as in Case 5 (Shearer et al., 1979) 
and for the same reasons. It should be mentioned also that the prediction of 
Cases 6 and 7 are obtained with the coefficients of the turbulence model given 
by Table 3.1. The optimized value for C^ in these two cases is equal to 2. 
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8.0 CONCLUSIONS AND RECOMMENDATIONS 

-The main objectives of the present study were as follows: 

1. to develop a mathematical model of turbulence for dilute two-phase flows 
starting from the exact transport equations of the turbulence kinetic 
energy and its dissipation rate; 

2. to develop a reliable formula for the calculation of the lateral 
diffusivity of heavy particles suspended in a homogeneous turbulent 
field, and 

3. to predict two-phase turbulent flows with phase changes based on modeled 
transport equations of mass, momentum of each phase, the concentration 
of vapor, and a two-equation turbulence model. 

In the presented model, the third-order correlations containing particle 
volume fraction fluctuations are retained. The numerical results for all the 
predicted cases showed that those third-order correlations are negligible 
compared with the second-order ones (two orders of magnitude less). This 
means that the present study has only one new empirical coefficient in its 
turbulence model (C £ ^). This coefficient is determined from one set of data 
(Case 1) and used very successfully in all other cases of the same 
experiment. A sensitivity study was conducted to investigate the influence of 
the value of C £ ^ on the model predictions. By changing the value of by 
10%, the maximum change in any radial profile is less than 3%. 

The study of the effects of the dispersed phase on the carrier phase flow 
properties, mean and fluctuating components, shows the following results: 
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1. The momentum interchange between the two phases which reflects the 
degree of disequilibrium between the phases, is a function of the 
dispersed phase properties such as droplet diameter, density, and mass 
loading ratio. In the case of heavy particles suspended in a turbulent 
gaseous media, the momentum interchange terms and all the corresponding 
turbulent correlations should be considered in the governing equations 
of both the mean motion and the turbulence model. 

2. The effect of partial or complete droplet evaporation is reflected on 
the velocity distribution of the different size groups. The smaller the 
diameter of the droplet is, the less the relative velocity between 
droplets and the fluid, and the higher the turbulent diffusivity of that 
group. 

3. Due to the co-existence of the dispersed phase and the carrier phase in 
the same control volume, a significant reduction in the turbulent shear 
stress and the kinetic energy of turbulence of the carrier phase is 
observed. The reduction in the turbulence energy or the increase in the 
dissipation rate of that energy is caused by the fluctuating relative 
velocity between the particles and the carrier phase and the turbulent 
correlation between this velocity and other fluctuating quantities such 
as volume fractions and carrier fluid velocity. The reduction in the 
kinetic energy of turbulence is proportional to the loading ratio but 
not linearly. 

A reliable expression for calculating the Schmidt number, defined as the 
ratio of particle diffusivity to fluid point diffusivity, of heavy particles 
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suspended in a turbulent flow is developed (Equations 5.30 to 5.33). The 

predictions using that _ formula are compared . with recent well-defined 

experimental data for the dispersion of a single particle. The agreement 
between the predictions and data is very good.’ 

Using the turbulence model presented in this work, predictions of the 
different cases for either solid particles or evaporating sprays, are 
generally in good agreement with the most recent well-defined experimental 
data. 

Further extension of the present work includes: 

1. obtaining optical measurements for the flow properties of the ideal 

spray experiment of Case 4 to validate the present model and to support 
the turbulent spray models in general. 


2. predicting a ducted recirculating turbulent two-phase flow (elliptic 
flow). The predictions should be compared with a well defined data set. 

3. predicting a ducted turbulent two-phase flow with heat transfer. The 
interaction between the evaporating droplet and the duct walls, and the 
heat and mass transfer to the wall, should be considered in the model. 
The density fluctuation effects should also be considered. 

4. predicting the dense portion of the spray. Droplet-droplet interaction 
effects, collision, and shattering must be considered in the model. 
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APPENDIX A 


Material Properties of the Spray 
Table A-l 

Physical Properties of Liquid Droplets* 


Property 

Methanol 

(ch 4 o) 

Freon .11 

(ccl 3 f) 

Saturated vapor pressure (P ), N/m 

0.207 


Latent heat of vaporization (L), KJ/Kg 

50.0 

181.32 

Density (p 2 >, Kg/m 3 

810 

1518 

Saturation temperature (T s ), °K 

292 

■: 240.3 

Boiling temperature (T^), °K 

347.71 

296.7 

Molecular weight, W v 

32 

137.37 

Viscosity of liquid material, 10^Kg/ms 

5.09 

4.05 

3 

Surface tension (y ) , 10 N/m 

21.8 

7.5 

Diffusivity of the evaporating material ( 5 ), 

10 5 m 2 /s 1.35 

2.85 


* Obtained at 30°C and P = 1 atm. (Vargaftik 1975 and ASHRAE 1969) 
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APPENDIX B " 


Modeled Transport Equations in Cartesian Tensor Notations 
Substituting in the time-averaged equations presented in the subsections 
2.3 & 3.3 by the modeling approximations for various turbulent correlations 
discussed in the subsection 3.4, the modeled transport equations in Cartesian 
tensor notations are obtained and will be presented in this Appendix. 


The continuity equations of the carrier phase 
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Third Order Correlation 
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The momentum equation of the k_ phase 
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The conentration equation 
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The turbulence kinetic energy equation (K) 
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APPENDIX C 


Initial Conditions of the Different Cases 
Table C-2. Experimental Flow Conditions of Modarress et al. (1984)* 


Gas-Phase (Air): 

Case 1 

Case 2 

Case 3 

Centerline velocity, U v (m/s) 

12.6 

12.6 

13.4 

Exponent, n, of power law velocity 
Profile U x /U x>c « (l-(2r/Q)) l /n 

C — 

6.6 — — 

> 

Turbulence Intensity 

< (0.04 

+ 0.1 r/D)— 

— — > 

Density, (Kg/m 3 ) 

< 

• 

00 

1 

1 

1 

1 


Mass flow rate m^ (Kg/s) 

3 .76x10” 3 

3.76xl0 -3 

4xl0 -3 

Reynolds number R e = (4m^/Try^D) 

13300 

13300 

14100 

Uniform mean velocity of surrounding 
stream, U x>s (m/s) 

< 

0.05 

> 

Intensity of turbulence in 
surrounding (u X(g / U x>s ) 

< 

0.1- 

> 


* Measured at 0.1D Downstream of Pipe Exit 
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Solid-Phase (Glass Beads): 

Case 1 

Case 2 

Case 3 

Particle diameter (microns) 

50 

50 

200 

Particle density, P 2 (Kg/m 3 ) 

< 

— 2990 

> 

Centerline velocity, V v _ (m/s) 

12.4 

12.4 

10.2 

Exponent, n, of power law velocity 
profile 

< 

27.6 

— > 

Mass flow rate iiL (Kg/s) 

1.2xl0 -3 

3.2xlO -3 

3 .2x1 0 - - 3 

Ratio of mass flow rates 

* = m„/m 
0 2 1 

0.32 

0.85 

0.8 

Ratio of volumetric fractions = 

V*1 ' <V i 'l )(D l U x,av. /c 2 V x,av. ) 

l.lxlO -4 

2.9xl0 -4 

3. 5 2x1 O' 4 
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Table 03: Experimental Flow Conditions of Shearer et al. (1979)* (Case 5) 
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Measured at 170 D Downstream of the Injection Nozzle 
S.ph. = single phase flow 



Table C-4. Experimental Flow Conditions of Solomon etal. (1984) 
Properties of the Carrier Phase at 
50 D Downstream of the Nozzle (Cases 6 & 7) 
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Table C-4b. Properties of the k cn Phases 
at 50 D Downstream of the Nozzle (Case 6) 
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Table C-l: The Considered Cases 


Case 

Number 

Dispensed 

Phase 

Material 

Diameter 

of 

Particles 
d ra 

Mass 

Loading 

Ratio 

X o 

Region 
- of 
Study 
z/D 

Reference 

1 

Glass 

50 

0.32 

0.1-20 

Modarress 
et al. (1984) 

2 

Glass 

50 

0.85 

0.1-20 

Modarress 
et al. (1984) 

3 

Glass 

200 

0.8 

Q. 1-20 

Modarress 
et al. (1984) 

4 

Methanol 

20-100 

0. 1-0.5 

0.1-20 

Idealized Flow 

5 

Freon- 1 1 

27 

6.88 

170-510 

. Shearer 
et al. (1979) 

6 

Freon- 1 1 

17.5-52.5 

7 .71 ; 

50-500 

Solomon 
et al. (1984) 

7 

Freon- l 1 

15-100 

15.78 

50-500 

Solomon 


. et al. (1984) 


Table C-4c. Mean Velocities of for k cn Phases 
at 50 D Downstream of the Nozzle (Case 7) 
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Table C-4d. Volume Fractions of the k cn Phases 
at 50 D Downstream of the Nozzle (Case 7) 
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